且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

用scipy拟合微分方程

更新时间:2022-06-26 23:59:08

您需要提供函数 f(t,b,c),该函数在中给出了一个参数或参数列表> t 返回参数处的函数值.这需要做一些工作,或者通过确定 t 的类型,或者通过使用一种可以以下两种方式起作用的构造:

You need to provide a function f(t,b,c) that given an argument or a list of arguments in t returns the value of the function at the argument(s). This requires some work, either by determining the type of t or by using a construct that works either way:

def f(t,b,c): 
    tspan = np.hstack([[0],np.hstack([t])])
    return scint.odeint(pend, y0, tspan, args=(b,c))[1:,0]

popt, pcov = scopt.curve_fit(f, t, sol[:,0], p0=guess)

返回 popt = array([0.25,5.]).

可以扩展它以容纳更多参数,

This can be extended to fit even more parameters,

def f(t, a0,a1, b,c): 
    tspan = np.hstack([[0],np.hstack([t])])
    return scint.odeint(pend, [a0,a1], tspan, args=(b,c))[1:,0]

popt, pcov = scopt.curve_fit(f, t, sol[:,0], p0=guess)

这会导致 popt = [3.04159267e + 00,-2.38543640e-07,2.49993362e-01,4.99998795e + 00] .

另一种可能性是显式计算目标解决方案的差异的平方范数,并将最小化应用于如此定义的标量函数.

Another possibility is to explicitly compute the square norm of the differences to the target solution and apply minimization to the so-defined scalar function.

 def f(param): 
     b,c = param
     t_sol = scint.odeint(pend, y0, t, args=(b,c))
     return np.linalg.norm(t_sol[:,0]-sol[:,0]);

res = scopt.minimize(f, np.array(guess))

res

      fun: 1.572327981969186e-08
 hess_inv: array([[ 0.00031325,  0.00033478],
                  [ 0.00033478,  0.00035841]])
      jac: array([ 0.06129361, -0.04859557])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 518
      nit: 27
     njev: 127
   status: 2
  success: False
        x: array([ 0.24999905,  4.99999884])