Подберите дифференциальное уравнение со Сципи - PullRequest
Новая
       11

Подберите дифференциальное уравнение со Сципи

0 голосов
/ 13 декабря 2018

как мне установить дифференциальную функцию из следующего учебника по scipy

Учебник по дифференциальному уравнению Scipy ?

В конце я хочу указать некоторые точки данных, следующие занабор из двух дифференциальных уравнений с шестью параметрами, но я хотел бы начать с простого примера.До сих пор я пробовал функции scipy.optimize.curve_fit и scipy.optimize.leastsq, но я никуда не попал.

Так вот, как далеко я зашел: 

import numpy as np
import scipy.optimize as scopt
import scipy.integrate as scint
import scipy.optimize as scopt

def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt


def test_pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

b = 0.25
c = 5.0

y0 = [np.pi - 0.1, 0.0]
guess = [0.5, 4]
t = np.linspace(0, 1, 11)

sol = scint.odeint(pend, y0, t, args=(b, c))

popt, pcov = scopt.curve_fit(test_pend, guess, t, sol)

со следующим сообщением об ошибке: 

ValueError: too many values to unpack (expected 2)

И я извиняюсь, потому что это довольно простой вопрос, но я не могу заставить его работать.Так что заранее спасибо.

1 Ответ

0 голосов
/ 13 декабря 2018

Вам необходимо предоставить функцию f(t,b,c), которая с аргументом или списком аргументов в t возвращает значение функции для аргумента (ов).Это требует некоторой работы, либо путем определения типа t, либо с помощью конструкции, которая работает в любом случае: 

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. ]).

Это может быть расширено доподходит еще больше параметров, 

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].

Другая возможность состоит в том, чтобы явно вычислить квадратную норму разностей для целевого решения и применить минимизацию копределенная скалярная функция. 

 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])
...