I would like to describe an experimental data ('exp=((10.0052, 87.9716), (11.0433, 86.3866), (12.359, 84.8015), (13.7127,83.2164),...'), where the first number in each pair is a temperature (T) and the second number is magnetic susceptibility (\chi) by the following formula:
$$a(T-T_C)\chi^{1/c}+b\chi^{1/c}*\chi^{1/d}H^{1/d}$$, where $a$,$b$,$c$,$d$,$T_C$ and $C$ are parameters (initial values of the parameters are $a=0.79$, $b=0.00893$, $c=1.38$, $d=0.428$, $T_C=315$, $C=1$) $T_C>0$, $C>0$
In my previous question @Reinderien did it for 4 parameters ($a$,$b$,$c$,$d$). The code is given below.
Could you show me please how to change the code so that the minimization is based on 6 parameters, i.e. add two more parameters ($TC$ and $C$) to the existing 4 parameters ($a$,$b$,$c$,$d$). Could you please plot a curve obtained using the parameters obtained during minimization together with the experimental data.
In the code C≡CC
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import NonlinearConstraint, minimize
exp=((10.0052, 87.9716), (11.0433, 86.3866), (12.359, 84.8015), (13.7127, 83.2164), (14.9939, 82.0276), (16.3367, 80.8388), (17.6282, 80.0463), (18.9464, 78.8575), (18.9934, 78.8575), (18.9969, 78.8575), (18.9962, 78.8575), (19.4221, 78.4612), (20.8034, 77.6687), (22.1306, 76.4799), (23.4352, 75.6873), (24.7503, 75.291), (26.0295, 74.4985), (27.3692, 73.706), (30.0264, 72.9134), (31.3262, 72.5172), (32.6313, 72.1209), (33.9404, 71.7246), (35.2494, 71.3284), (36.551, 70.9321), (37.8445, 70.5358), (40.5813, 69.7433), (41.7696, 69.347), (43.0624, 69.347), (44.3609, 68.9507), (45.6536, 68.5545), (46.9616, 68.5545), (48.254, 68.1582), (49.5428, 67.7619), (52.1298, 67.3657), (53.4147, 67.3657), (54.7091, 66.9694), (56.0001, 66.9694), (57.2883, 66.5731), (58.5865, 66.5731), (59.8911, 66.1769), (62.4856, 66.1769), (63.7665, 65.7806), (65.0689, 65.7806), (66.3588, 65.3843), (67.6581, 65.3843), (68.9546, 65.3843), (70.2572, 64.9881), (72.8482, 64.9881), (74.1637, 64.9881), (75.4546, 64.5918), (76.7495, 64.5918), (78.0347, 64.5918), (79.3376, 64.1955), (80.6437, 64.1955), (83.5861, 64.1955), (84.8957, 63.7993), (86.2123, 63.7993), (87.5317, 63.7993), (88.8468, 63.7993), (90.1674, 63.7993), (91.4775, 63.403), (94.2769, 63.403), (95.5863, 63.403), (96.9007, 63.403), (98.2013, 63.403), (99.523, 63.0067), (100.839, 63.0067), (102.147, 63.0067), (104.812, 63.0067), (106.131, 63.0067), (107.448, 62.6104), (108.789, 62.6104), (110.099, 62.6104), (111.414, 62.6104), (112.756, 62.6104), (115.404, 62.6104), (116.737, 62.6104), (118.059, 62.6104), (119.376, 62.2142), (120.7, 62.2142), (122.025, 62.2142), (123.367, 62.2142), (126.006, 62.2142), (127.352, 62.2142), (128.674, 62.2142), (129.999, 62.2142), (131.333, 62.2142), (132.652, 61.8179), (133.987, 61.8179), (136.652, 61.8179), (137.987, 61.8179), (139.324, 61.8179), (140.645, 61.8179), (141.989, 61.8179), (143.316, 61.8179), (144.646, 61.8179), (147.323, 61.8179), (148.651, 61.4216), (149.978, 61.4216), (151.303, 61.4216), (152.637, 61.4216), (153.949, 61.4216), (155.268, 61.4216), (157.929, 61.4216), (159.252, 61.4216), (160.581, 61.0254), (161.889, 61.0254), (163.219, 61.0254), (164.546, 61.0254), (165.873, 61.0254), (168.524, 61.0254), (169.852, 60.6291), (171.173, 60.6291), (172.494, 60.6291), (173.83, 60.6291), (175.16, 60.6291), (176.473, 60.6291), (179.122, 60.2328), (180.459, 60.2328), (181.78, 60.2328), (183.106, 59.8366), (184.439, 59.8366), (185.757, 59.8366), (187.079, 59.8366), (189.73, 59.4403), (191.06, 59.4403), (192.386, 59.044), (193.717, 59.044), (195.037, 59.044), (196.372, 58.6478), (197.693, 58.6478), (200.361, 58.2515), (201.689, 57.8552), (203.008, 57.8552), (204.335, 57.459), (205.651, 57.459), (206.996, 57.0627), (208.324, 57.0627), (210.977, 56.6664), (212.292, 56.2702), (213.602, 55.8739), (214.947, 55.8739), (216.259, 55.4776), (217.583, 55.0813), (218.895, 55.0813), (221.56, 54.2888), (222.865, 53.8925), (224.195, 53.4963), (225.505, 53.1), (226.824, 53.1), (228.192, 52.7037), (229.508, 52.3075), (232.191, 51.5149), (233.507, 51.1187), (234.836, 50.7224), (236.161, 50.3261), (237.479, 49.9299), (238.807, 49.5336), (240.13, 48.741), (242.774, 47.9485), (244.096, 47.5522), (245.417, 47.156), (246.74, 46.7597), (248.057, 46.3634), (249.382, 45.5709), (250.457, 46.2508), (251.117, 46.0365), (252.897, 45.4166), (253.112, 45.3033), (254.115, 44.9356), (255.112, 44.5707), (256.107, 44.2179), (257.936, 43.554), (259.077, 43.1208), (260.395, 42.633), (260.555, 42.5688), (261.881, 42.0749), (263.515, 41.4803), (264.389, 41.108), (264.559, 41.047), (264.724, 40.9939), (264.891, 40.9305), (265.058, 40.8548), (265.219, 40.807), (265.384, 40.7401), (265.549, 40.6799), (265.707, 40.6116), (265.872, 40.5381), (266.041, 40.4919), (266.218, 40.4131), (266.379, 40.365), (266.532, 40.2826), (266.705, 40.2386), (266.876, 40.1762), (267.038, 40.1109), (267.201, 40.0532), (267.373, 39.9972), (267.543, 39.9247), (267.71, 39.8555), (267.874, 39.7915), (268.035, 39.7336), (268.194, 39.6811), (268.367, 39.6056), (268.542, 39.5387), (268.705, 39.4762), (268.866, 39.3997), (269.029, 39.3423), (269.201, 39.2861), (269.374, 39.22), (269.539, 39.1538), (269.698, 39.0823), (269.87, 39.0231), (270.036, 38.9709), (270.194, 38.8925), (270.371, 38.8394), (270.536, 38.7638), (270.697, 38.7051), (270.868, 38.6425), (271.038, 38.574), (271.208, 38.5115), (271.367, 38.4455), (271.54, 38.3782), (271.711, 38.3156), (271.868, 38.2524), (272.03, 38.1917), (272.2, 38.1262), (272.361, 38.0706), (272.525, 38.0021), (274.156, 37.3821), (274.229, 37.3188), (274.362, 37.27), (274.524, 37.2048), (274.705, 37.1411), (274.877, 37.0827), (275.029, 37.0104), (275.196, 36.9607), (275.37, 36.8851), (275.526, 36.8172), (275.683, 36.7629), (275.855, 36.6982), (276.025, 36.6346), (276.185, 36.5695), (276.345, 36.5063), (276.506, 36.4347), (276.68, 36.3789), (276.845, 36.315), (277.012, 36.2632), (277.184, 36.1689), (277.35, 36.1295), (277.517, 36.0462), (277.676, 35.9975), (277.839, 35.9304), (278.009, 35.8579), (278.18, 35.8111), (278.341, 35.7453), (278.5, 35.6707), (278.67, 35.6269), (278.836, 35.5542), (279.005, 35.486), (279.176, 35.4266), (279.333, 35.3622), (279.502, 35.2938), (279.671, 35.2304), (279.833, 35.1688), (279.997, 35.1014), (280.169, 35.0437), (280.337, 34.9787), (280.498, 34.9181), (280.67, 34.8453), (280.837, 34.7763), (280.995, 34.7195), (281.164, 34.6492), (281.332, 34.5889), (281.5, 34.523), (281.664, 34.4628), (281.83, 34.3953), (282.004, 34.3438), (282.164, 34.2835), (282.324, 34.214), (282.493, 34.1503), (282.67, 34.087), (282.833, 34.0348), (282.996, 33.9549), (283.164, 33.9051), (284.809, 33.291), (284.864, 33.2491), (284.993, 33.1839), (285.16, 33.111), (285.323, 33.0541), (285.487, 32.9876), (285.655, 32.9227), (285.821, 32.8639), (285.991, 32.8064), (286.164, 32.7342), (286.329, 32.6806), (286.494, 32.606), (286.66, 32.562), (286.824, 32.4928), (286.985, 32.4283), (287.149, 32.3634), (287.321, 32.2873), (287.487, 32.2322), (287.644, 32.1765), (287.817, 32.1208), (287.988, 32.0552), (288.146, 31.9866), (288.306, 31.9345), (288.476, 31.8621), (288.646, 31.8045), (288.815, 31.7536), (288.985, 31.6971), (289.14, 31.6201), (289.307, 31.551), (289.481, 31.4893), (289.645, 31.444), (289.81, 31.3815), (289.975, 31.311), (290.141, 31.252), (290.306, 31.1829), (290.466, 31.1326), (290.636, 31.0702), (290.807, 31.011), (290.967, 30.9446), (291.139, 30.8823), (291.304, 30.8283), (291.47, 30.7563), (291.646, 30.706), (291.809, 30.6516), (291.974, 30.578), (292.141, 30.5155), (292.306, 30.4571), (292.47, 30.3986), (292.634, 30.3283), (292.805, 30.2715), (292.971, 30.212), (293.127, 30.1523), (293.299, 30.0891), (293.475, 30.0263), (293.634, 29.9634), (293.794, 29.9016), (295.41, 29.3469), (295.467, 29.2859), (295.612, 29.2319), (295.787, 29.1886), (295.954, 29.1198), (296.116, 29.0857), (296.277, 29.0178), (296.445, 28.9598), (296.617, 28.9062), (296.78, 28.8415), (296.946, 28.7719), (297.111, 28.7037), (297.277, 28.6442), (297.444, 28.5661), (297.603, 28.5065), (297.772, 28.438), (297.945, 28.3781), (298.11, 28.3183), (298.27, 28.2553), (298.436, 28.1991), (298.608, 28.1362), (298.775, 28.089), (298.934, 28.0215), (299.089, 27.9536), (299.26, 27.8898), (299.431, 27.8379), (299.596, 27.7736), (299.77, 27.7224), (299.939, 27.6707), (300.1, 27.6075), (300.268, 27.5584), (300.436, 27.4834), (300.595, 27.445), (300.761, 27.3805), (300.928, 27.3191), (301.085, 27.2599), (301.248, 27.2049), (301.416, 27.1372), (301.582, 27.0904), (301.742, 27.0381), (301.902, 26.9595), (302.077, 26.9182), (302.254, 26.8639), (302.409, 26.7965), (302.573, 26.751), (302.744, 26.6902), (302.904, 26.6289), (303.075, 26.5635), (303.24, 26.5116), (303.398, 26.4599), (303.571, 26.3952), (303.75, 26.3447), (303.915, 26.2879), (304.07, 26.2416), (304.236, 26.1701), (304.396, 26.1088))
def CH(
chi: np.ndarray, T: float, H: float, TC: float, CC: float, a: float, b: float, c: float, d: float,
) -> np.ndarray:
err = (
a*(T - TC)*chi**(1/c)
+ b * chi**(1/c) * chi**(1/d) * H**(1/d)
+ CC/(T-TC)-chi
- 1
)
return err
exp_T, exp_chi = np.array(exp).T
def dchi(params: np.ndarray) -> float:
chi = params[2:]
chi_err = chi - exp_chi
err = chi_err.dot(chi_err) #squre errors
return err
def dchi_jac(params: np.ndarray) -> np.ndarray:
chi = params[2:]
jac = np.zeros_like(params)
jac[2:] = 2*(chi - exp_chi)
return jac
def dchi_root(params: np.ndarray) -> np.ndarray:
# for each value of exp_T, there is some unknown optimal value of chi such that
# CH(chi, T) == 0
abcdCTC, chi = np.split(params, (6,))
a, b, c, d ,CC, TC = abcdCTC
root_err = CH(
chi=chi, T=exp_T, a=a, b=b, c=c, d=d, TC=320, CC=1, H=1000,
)
return root_err
exp_T, exp_chi = np.array(exp).T
first_guess = np.concatenate((
(0.79, 0.00893, 0.428, 1.38, 320, 1),
exp_chi,
))
result = minimize(
fun=dchi, jac=dchi_jac,
x0=first_guess,
constraints=NonlinearConstraint(
fun=dchi_root,
lb=0, ub=0,
),
options={'maxiter': 10, 'disp': True},
)
print(result.message)
print(f'Chi error: {result.fun:.3f}')
print(f'Root error: {np.abs(dchi_root(result.x)).sum():.3f}')
print(f'Parameters: {result.x[:6]}')
chi_optimized = result.x[2:]
plt.figure(figsize=(10, 5))
plt.plot(exp_T, exp_chi, 'o', label='Experimental chi')
plt.xlabel('T')
plt.ylabel('chi')
plt.legend()
plt.grid(True)
plt.title('Graph of chi function')
plt.show()
the error
Traceback (most recent call last):
File "C:\Users\...\PycharmProjects\pythonProject7\main.py", line 50, in <module>
result = minimize(
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_minimize.py", line 719, in minimize
res = _minimize_slsqp(fun, x0, args, jac, bounds,
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_slsqp_py.py", line 374, in _minimize_slsqp
sf = _prepare_scalar_function(func, x, jac=jac, args=args, epsilon=eps,
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_optimize.py", line 383, in _prepare_scalar_function
sf = ScalarFunction(fun, x0, args, grad, hess,
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 158, in __init__
self._update_fun()
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 251, in _update_fun
self._update_fun_impl()
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 155, in update_fun
self.f = fun_wrapped(self.x)
File "C:\Users\...\PycharmProjects\pythonProject7\venv\lib\site-packages\scipy\optimize\_differentiable_functions.py", line 137, in fun_wrapped
fx = fun(np.copy(x), *args)
File "C:\Users\...\PycharmProjects\pythonProject7\main.py", line 24, in dchi
chi_err = chi - exp_chi
ValueError: operands could not be broadcast together with shapes (400,) (396,)
You were using a couple of incorrect indices, and failed to pass the value for
TCandCC(which should just be namedC). You also have incorrect notation forcandd, which have different variable names in your cited paper - but I have not modified this.