DjangoBB LoFi version

Полная версия: Нелинейный МНК с интегралами

Начало » Центр помощи » Нелинейный МНК с интегралами

909404

Авг. 22, 2017 22:33:08

Всем привет!
ВОзникли трудности с подбором оптимального значения параметра на питоне, гугл не помог, пишу сюда.

Есть массив точек Cv - T.

ССылка на картинку…

На картинке показана модель Дебая для твердого тела. N, k- константы. Td - подбираемый параметр. Соотвественно, имея достаточное количество пар точек Cv - T можно подобрать оптимальное значение параметра Td.

Код выдает ошибку
ValueError: operands could not be broadcast together with shapes (2,) (75,)

На сколько я понимаю, в функции D(T,Td) T - список экспериментальных точек, а Td - число.
Вполне вероятно, что я использую не тот функционал. Буду благодарен за любую помощь!

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy import optimize

#Ввод экспериментальных данных
exp_data = np.loadtxt('Cp.txt')

#Constants used in calculations
h = 1.05E-34 
k = 1.38E-23 
Na = 6.02E+23

T_exp = exp_data[:,0] #x-values
Cp_exp = exp_data[:,1] #y-values

def Cp(T,Td,a0):
    return a0*T*(6*Na*k*D(T, Td))**2 + 6*Na*k*D(T, Td)

def integrand(x):
    return (x**4)*np.exp(x)/((np.exp(x)-1)**2)

def D(T, Td):
    y, err = integrate.quad(integrand, 0.0, Td/T)
    return 3*(T/Td)**3*y

def residual(Td, a0, T):
    return Cp_exp - Cp(T, Td, a0)

p0 = [2E+13, -2]
popt, pcov = optimize.leastsq(residual, p0, args=(T_exp, Cp_exp)

scidam

Авг. 23, 2017 04:48:32

Поскольку я не специалист в статистической физике, до конца я задачу не понял: т.к. исходная формула выглядит несколько по-другому, чем код. В коде есть h, a0 да и Cp Cv путают меня… Но, думаю следующий работающий код, поможет вам достичь того, чего вы хотите:

 import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from scipy import optimize
#Ввод экспериментальных данных
# exp_data = np.loadtxt('Cp.txt')
exp_data = 3 + np.random.rand(75, 2) # exp. data simulation
#Constants used in calculations
h = 1.05E-34 # unused var.
k = 1.38E-23
Na = 6.02E+23
T_exp = exp_data[:,0] #x-values
Cp_exp = exp_data[:,1] #y-values
# def Cp(T,Td,a0):
#     return a0*T*(6*Na*k*D(T, Td))**2 + 6*Na*k*D(T, Td)
def Cp(T, Td): #Given formula doesn't depend on <a0>
     return Na * k * D(T, Td)
def integrand(x): # be pep8 compliant...
    return (x ** 4) * np.exp(x) / ((np.exp(x) - 1) ** 2)
@np.vectorize  # T is an array, but Td is a scalar value; we need to vectorize D!!! Important!
def D(T, Td): # be pep8 compliant...
    y, err = integrate.quad(integrand, 1.0e-10, Td / T)
    return 9 * y * (T / Td) ** 3
def residuals(Td): # the only one parameter Td is estimated!?
    return Cp_exp - Cp(T_exp, Td)
initital_Td = 1.0
popt, pcov = optimize.leastsq(residuals, initital_Td)
print(popt, pcov)

909404

Авг. 25, 2017 20:18:45

Спасибо, добрый человек!
Все заработало, и я сильно продвинулся вперед!