我正在努力改变我的代码来解决一个常微分方程系统,从初值问题到边值问题。我已经试过很多次了,但我认为我犯的错误在逻辑上是不正确的。因此,我没有粘贴更改代码,而是将其粘贴到我的原始代码下面,这很好
下面提到的代码用于解决一个具有函数的ODE系统,然后我使用粒子群优化(PSO)算法进行优化。我想使用函数的相同方程,用边界条件t(0)=1和t(1)=2求解_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)
如果你有一首像这样的颂歌
你可以把它当作
或作为
通过编码必要的函数和初始猜测,可以切换到边界问题
请注意,如果您的新解算器函数版本支持以与
odeint
相同的方式传递参数,则必须进行尝试。如果不可能,请使用lambda构造作为包装器或显式定义的函数相关问题 更多 >
编程相关推荐