Как минимизировать функцию по нескольким параметрам? (часть 2) ⇐ Python
-
Anonymous
Как минимизировать функцию по нескольким параметрам? (часть 2)
Я хотел бы описать экспериментальные данные ('exp=((10.0052, 87.9716), (11.0433, 86.3866), (12.359, 84.8015), (13.7127,83.2164),...'), где первое число в каждой паре стоит температура (Т), а второе число — магнитная восприимчивость (\chi) по следующей формуле:
$$a(T-T_C)\chi^{1/c}+b\chi^{1/c}*\chi^{1/d}H^{1/d}$$, где $ a$,$b$,$c$,$d$,$T_C$ и $C$ — параметры (начальные значения параметров $a=0,79$, $b=0,00893$, $c=1,38$, $ d=0.428$, $T_C=315$, $C=1$) $T_C>0$, $C>0$
В моем предыдущем вопросе @Reinderien сделал это для 4 параметров ($a$,$b$,$c$,$d$). Код приведен ниже.
Не могли бы вы показать мне, как изменить код, чтобы минимизация основывалась на 6 параметрах, т.е. к существующим 4 параметрам ($a$,$b$) добавить еще два параметра ($TC$ и $C$) ,$c$,$d$). Не могли бы вы построить кривую, полученную с использованием параметров, полученных при минимизации, вместе с экспериментальными данными.
импортировать numpy как np импортировать matplotlib.pyplot как plt из scipy.optimize импортировать NonlinearConstraint, минимизировать 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.32) 84), (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.9) 881), (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.2 123, 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.6) 104), (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.36) 7, 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.4) 216), (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.83) 66), (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.2) 702), (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.49) 63), (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.570) 9), (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.5) 54), (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.9) 247), (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.83) 94), (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.315) 6), (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, 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[4:] chi_err = chi - exp_chi err = chi_err.dot(chi_err) return err def dchi_jac(params: np.ndarray) -> np.ndarray: chi = params[4:] jac = np.zeros_like(params) jac[4:] = 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 abcd, chi = np.split(params, (4,)) a, b, c, d = abcd root_err = CH( chi=chi, T=exp_T, a=a, b=b, c=c, d=d, H=1000, TC=315, C=2, ) return root_err exp_T, exp_chi = np.array(exp).T first_guess = np.concatenate(( (0.79, 0.00893, 0.428, 1.38), 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[:4]}') 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()
Я хотел бы описать экспериментальные данные ('exp=((10.0052, 87.9716), (11.0433, 86.3866), (12.359, 84.8015), (13.7127,83.2164),...'), где первое число в каждой паре стоит температура (Т), а второе число — магнитная восприимчивость (\chi) по следующей формуле:
$$a(T-T_C)\chi^{1/c}+b\chi^{1/c}*\chi^{1/d}H^{1/d}$$, где $ a$,$b$,$c$,$d$,$T_C$ и $C$ — параметры (начальные значения параметров $a=0,79$, $b=0,00893$, $c=1,38$, $ d=0.428$, $T_C=315$, $C=1$) $T_C>0$, $C>0$
В моем предыдущем вопросе @Reinderien сделал это для 4 параметров ($a$,$b$,$c$,$d$). Код приведен ниже.
Не могли бы вы показать мне, как изменить код, чтобы минимизация основывалась на 6 параметрах, т.е. к существующим 4 параметрам ($a$,$b$) добавить еще два параметра ($TC$ и $C$) ,$c$,$d$). Не могли бы вы построить кривую, полученную с использованием параметров, полученных при минимизации, вместе с экспериментальными данными.
импортировать numpy как np импортировать matplotlib.pyplot как plt из scipy.optimize импортировать NonlinearConstraint, минимизировать 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.32) 84), (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.9) 881), (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.2 123, 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.6) 104), (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.36) 7, 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.4) 216), (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.83) 66), (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.2) 702), (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.49) 63), (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.570) 9), (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.5) 54), (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.9) 247), (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.83) 94), (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.315) 6), (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, 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[4:] chi_err = chi - exp_chi err = chi_err.dot(chi_err) return err def dchi_jac(params: np.ndarray) -> np.ndarray: chi = params[4:] jac = np.zeros_like(params) jac[4:] = 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 abcd, chi = np.split(params, (4,)) a, b, c, d = abcd root_err = CH( chi=chi, T=exp_T, a=a, b=b, c=c, d=d, H=1000, TC=315, C=2, ) return root_err exp_T, exp_chi = np.array(exp).T first_guess = np.concatenate(( (0.79, 0.00893, 0.428, 1.38), 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[:4]}') 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()
