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 TC and CC (which should just be named C). You also have incorrect notation for c and d, which have different variable names in your cited paper - but I have not modified this.
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))
exp_T, exp_chi = np.array(exp).T
def CH(
chi: np.ndarray, T: float, H: float, TC: float, C: 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)
+ C/(T - TC)
- chi
- 1
)
return err
def dchi(params: np.ndarray) -> float:
chi = params[6:]
chi_err = chi - exp_chi
err = chi_err.dot(chi_err) # least-squared error
return err
def dchi_jac(params: np.ndarray) -> np.ndarray:
chi = params[6:]
jac = np.zeros_like(params)
jac[6:] = 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
(a, b, c, d, C, TC), chi = np.split(params, (6,))
root_err = CH(
chi=chi, T=exp_T, a=a, b=b, c=c, d=d, TC=TC, C=C, H=1000,
)
return root_err
def solve() -> None:
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,
},
)
params, chi_optimized = np.split(result.x, (6,))
print(result.message)
print(f'Chi error: {result.fun:.3f}')
print(f'Root error: {np.abs(dchi_root(result.x)).sum():.3f}')
print('Parameters:', params)
fig, ax = plt.subplots()
ax.plot(exp_T, exp_chi, label='Experiment')
ax.set_xlabel('T')
ax.set_ylabel('chi')
ax.legend()
ax.grid(True)
plt.show()
if __name__ == '__main__':
solve()
Iteration limit reached (Exit mode 9)
Current function value: 3354.871163641425
Iterations: 100
Function evaluations: 199
Gradient evaluations: 100
Iteration limit reached
Chi error: 3354.871
Root error: 0.001
Parameters: [ 6.10888764e-03 3.67890944e-05 1.88619221e+00 9.27603979e-01
3.52451910e+02 -1.11299938e+02]