将代码更改为ODE Python的边值问题

问题描述 投票:0回答:1

我正在努力更改代码以解决(常微分方程)ODE系统从初值问题到边值问题的问题。我已经尝试了很多次,但是我认为我犯的错误在逻辑上是不正确的。因此,我粘贴了正常的原始代码,而不是粘贴更改代码。

下面提到的代码用于解决具有功能odeint]的ODE系统,然后使用粒子群优化(PSO)算法进行优化。我想在函数边界条件为t(0)= 1和t(1)= 2的情况下使用函数solve_bvp的相同方程式。下面提到的是代码。谢谢

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)

我正在努力更改代码以解决(常微分方程)ODE系统从初值问题到边值问题的问题。我已经尝试了很多次,但我想我正在做...

python numpy scipy ode odeint
1个回答
0
投票

如果您有类似的ODE,则>

def f_ode(t,u): return [ u[1], -u[0] ]
© www.soinside.com 2019 - 2024. All rights reserved.