Изменение кода на граничную задачу для ODE Python - PullRequest
Новая
       85

Изменение кода на граничную задачу для ODE Python

0 голосов
/ 03 февраля 2020

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

Приведенный ниже код используется для решения системы ODE с функцией odeint , а затем я использую алгоритм Particle Swarm Optimization (PSO) для процесса оптимизации. Я хочу использовать те же уравнения с функцией solve_bvp с граничными условиями t (0) = 1 и t (1) = 2 , Ниже упоминается код. Спасибо 

from scipy import *
from scipy.integrate import odeint
from operator import itemgetter
import matplotlib
matplotlib.use('Agg')
from matplotlib.ticker import FormatStrFormatter
from pylab import *
from itertools import product
import itertools
from numpy import zeros_like
import operator
from pyswarm import pso

modelsOne = []
modelsTwo = []
modelsThree = []

#step 1 start## To build model structure library.  HIV model is a three-variable model, so we need three model structure liararys: modelsOne, modelsTwo, modelsThree.
              #  the model structure library contains all possible structures of the model to be sought.
def ModelsProduct(modelsOne, modelsTwo, modelsThree):
    modelsStepOne = list(product("+-",repeat = 4))
    modelsStepThree = [('a','a'),('a','b'),('a','c'),('b','b'),('b','c'),('c','c')]

    #produce modelsOne
    modelsStepTwo = [('b',),('c',)]
    for one in modelsStepOne:
        for two in modelsStepTwo:
            for three in modelsStepThree:
                modelsOne.append(one+two+three)

    #produce modelsTwo
    modelsStepTwo = [('a',),('c',)]
    for one in modelsStepOne:
        for two in modelsStepTwo:
            for three in modelsStepThree:
                modelsTwo.append(one+two+three)

    #produce modelsThree
    modelsStepTwo = [('a',),('b',)]
    for one in modelsStepOne:
        for two in modelsStepTwo:
            for three in modelsStepThree:
                modelsThree.append(one+two+three)
    return modelsOne, modelsTwo, modelsThree

modelsOne, modelsTwo,modelsThree = ModelsProduct(modelsOne, modelsTwo, modelsThree)



#step 1 end##



VarList = ["a","b","c"]
initial_condi = [100, 150, 50000]

dictVar = {'a':0, 'b': 1, 'c': 2}
ops = { "+": operator.add, "-": operator.sub }
t_range = arange(0.0,60.0,1.0)


def odeFunc(Y, t, x,dictVar):
    if x[-3] == 192:
        temp1 = 191
    else:
        temp1 = int(x[-3])
    if x[-2] == 192:
        temp2 = 191
    else:
        temp2 = int(x[-2])
    if x[-1] == 192:
        temp3 = 191
    else:
        temp3 = int(x[-1])
    modelOne = modelsOne[temp1]
    modelTwo = modelsTwo[temp2]
    modelThree = modelsThree[temp3]
    return GenModel(Y, x, modelOne,modelTwo,modelThree, dictVar)


def GenModel(Y,x,modelOne,modelTwo,modelThree, dictVar):
    dydt = zeros_like(Y)
    dydt[0] = ops[modelOne[0]](dydt[0],x[0]*Y[0])
    dydt[0] = ops[modelOne[1]](dydt[0],x[1]*Y[dictVar[modelOne[-3]]])
    dydt[0] = ops[modelOne[2]](dydt[0],x[2]*Y[dictVar[modelOne[-2]]]*Y[dictVar[modelOne[-1]]])
    dydt[0] = ops[modelOne[3]](dydt[0],x[3])

    dydt[1] = ops[modelTwo[0]](dydt[1],x[4]*Y[1])
    dydt[1] = ops[modelTwo[1]](dydt[1],x[5]*Y[dictVar[modelTwo[-3]]])
    dydt[1] = ops[modelTwo[2]](dydt[1],x[6]*Y[dictVar[modelTwo[-2]]]*Y[dictVar[modelTwo[-1]]])
    dydt[1] = ops[modelTwo[3]](dydt[1],x[7])

    dydt[2] = ops[modelThree[0]](dydt[2],x[8]*Y[2])
    dydt[2] = ops[modelThree[1]](dydt[2],x[9]*Y[dictVar[modelThree[-3]]])
    dydt[2] = ops[modelThree[2]](dydt[2],x[10]*Y[dictVar[modelThree[-2]]]*Y[dictVar[modelThree[-1]]])
    dydt[2] = ops[modelThree[3]](dydt[2],x[11])

    return dydt

## equations
def pendulum_equations(w, t):
    T, I, V = w
    dT = 80 - 0.15*T - 0.00002*T*V
    dI = 0.00002*T*V - 0.55*I
    dV = 900*0.55*I - 5.5*V - 0.00002*T*V
    return  dT, dI, dV

result_init = odeint(pendulum_equations, initial_condi, t_range)
result_init[:,2] =  result_init[:,2]/100

def myfunc(xRand):
    result_new = odeint(odeFunc, initial_condi, t_range, args=(xRand,dictVar))
    result_new[:,2] = result_new[:,2]/100
    result_sub = result_new - result_init
    return sum(result_sub*result_sub)


x = (0.15,0,0.00002,80,0.55,0,0.00002,0,5.5,495,0.00002,0,122,98,128)


lb = [0]*15
ub = [1,1,0.5,200,1,1,0.5,200,10,1000,0.5,200,192,192,192]


xopt1, fopt1 = pso(myfunc, lb, ub,omega= 0.7298,phip=1.49618,phig=1.49618,maxiter=1000,swarmsize= 1000,minstep=1e-20,minfunc=1e-20,debug = True)

1 Ответ

0 голосов
/ 03 февраля 2020

Если у вас есть ODE типа 

def f_ode(t,u): return [ u[1], -u[0] ]

, который вы можете решить как 

tspan = np.linspace(0,1,51);
u_init = [1.0, 0.0]
u = odeint(f_ode, u_init, tspan, tfirst=True)

или как 

res = solve_ivp(f_ode, tspan[[0,-1]], u_init, t_eval=tspan)
if res.success:
    u=res.y

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

def f_bc(u0, u1): return [u0[0]-1, u1[0]-2]
t = np.linspace(0,1,11);
u = [ 1+t, 1+0*t]
res = solve_bvp(f_ode,f_bc,t,u)
if res.success:
    u = res.sol(tspan)

Обратите внимание, что вы должны проверить, поддерживает ли ваша версия новых решающих функций передачу параметров так же, как и odeint. Если это невозможно, используйте лямбда-конструкции в качестве оболочек или явно определенные функции.

...