Одновременное численное соответствие двух уравнений методом наименьших квадратов - PullRequest
Новая
       74

Одновременное численное соответствие двух уравнений методом наименьших квадратов

0 голосов
/ 30 сентября 2018

Я пытаюсь привести в соответствие два упомянутых ниже уравнения, используя метод python leastsq, но не уверен, что это правильный подход.Первое уравнение содержит неполную гамма-функцию, а второе - немного сложное, и наряду с экспоненциальной функцией содержит член, который получается с помощью отдельной формулы подбора. 

J_mg = T_incomplete(hw/T_mag)
J_nmg = e^(-hw/T)*g(w,T)

Здесь g является функциейw и T и рассчитывается по заданной формуле соответствия.Я выполняю шаги, описанные в этом вопросе.

Вот что я сделал 

import numpy as np
from scipy.optimize import leastsq
from scipy.special import gammaincc
from scipy.special import gamma
from matplotlib.pyplot import plot


# generating data

NPTS = 10
hw = np.linspace(0.5, 10, NPTS)
j1 = np.linspace(0.001,10,NPTS)
j2 = np.linspace(0.003,10,NPTS)
T_mag = np.linspace(0.3,0.5,NPTS)

#defining functions

def calc_gaunt_factor(hw,T):
    fitting_coeff= np.loadtxt('fitting_coeff.txt', skiprows=1)
    #T is in KeV
    #K_b = 8.6173303(50)e−5 ev/K

    g = 0
    gamma = 0.0136/T
    theta= hw/T
    A= (np.log10(gamma**2) +0.5)*0.4
    B= (np.log10(theta)+1.5)*0.4
    for i in range(11):
        for j in range(11):
            g_ij = fitting_coeff[i][j]*(A**i)*(B**j)
            g = g_ij+g        
    return g

def j_w_mag(hw,T_mag):
    order= 0.001
    return np.sqrt(1/T_mag)*gamma(order)*gammaincc(order,hw/T_mag)

def j_w_nonmag(hw,T):
    gamma = 0.0136/T
    theta= hw/T
    return np.sqrt(1/T)*np.exp((-hw)/T)*calc_gaunt_factor(hw,T)

def residual_func(T,T_mag,hw,j1,j2):
    err_unmag = np.nan_to_num(j1 - j_w_nonmag(hw,T))
    err_mag = np.nan_to_num(j2 - j_w_mag(hw,T_mag))
    err= np.concatenate((err_unmag, err_mag))
    return err


par_init = np.array([.35])

best, cov, info, message, ler = leastsq(residual_func,par_init,args=(T_mag,hw,j1,j2),full_output=True)

print("Best-Fit Parameters:")
print("T=%s" %(best[0]))

Я получаю странное значение для моего подходящего параметра, Т.это правильный подход?Спасибо.

...