常微分⽅程求解器ODEsolver
另一种温暖⽂章⽬录
IVP 问题
对于常微分⽅程的初值问题
有以下解的存在性、唯⼀性定理：
假设定义在集合并且初值，函数对于变量是Lipschitz连续，则在与之间存在，使得初值问题
有唯⼀解。⽽且，如果在上是Lipschitz连续，则其在上存在唯⼀解。定义⼀个ODEsolver类来使⽤各种⽅法求解常微分⽅程，使⽤待求解的函数和初始值来初始化该类：class ODEsolver :
大硬盘def __init__(lf ,fun ,t0,y0):
lf .fun =fun
lf .t0=t0
lf .y0=y0
Euler 显式格式
Euler⽅法为：
代码为：
def eulersolver (lf ,tend ,step =None ):
if step is None :
step =(tend -lf .t0)/1e3
from ArrayTable import ArrayTable
result =ArrayTable (2,0)
result .tTableHeader (["$t$","$y(t)$"])
result .tTableUnit (["/","/"])
努力奋斗
result .append ([lf .t0,lf .y0])
t =lf .t0
while t <tend :
ynext =result .table [1].data [-1]+lf .fun (t ,result .table [1].data [-1])*step
result .append ([t +step ,ynext ])
台湾的风景名胜t +=step
return result
{=f (t ,y ),t ∈a ,b dt dy []y a =y ()0(1)
f (t ,y )[a ,b ]×[α,β]α<y <0βy a b c ⎩⎨⎧=f (t ,y ),t ∈a ,b dt dy []y a =y ()0t ∈[a ,c ](2)
y (t )f [a ,b ]×[−∞,+∞][a ,b ]{w =y 00w =w +hf t ,w i +1i (i i )(2)
隐式Euler 法
隐式Euler法的原理为：
迭代开始前⾸先要预估⼀个的初始值，⼀般有两种选择，使⽤作为初始值或者⽤显⽰⽅法先预估⼀个值，这⾥⽤预报矫正法来
估计其值。
def implicitEuler (lf ,tend , step =None ):
if step is None :
step = (tend - lf .t0) / 1e3
from ArrayTable import ArrayTable
result = ArrayTable (2, 0)
result .tTableHeader (["$t$", "$y(t)$"])
result .tTableUnit (["/", "/"])
result .append ([lf .t0, lf .y0])
t = lf .t0
while t <tend :
temp =lf .fun (t ,result .table [1].data [-1])
ynextinit =result .table [1].data [-1]+step /2.*(temp +lf .fun (t +step ,result .table [1].data [-1] + temp * step )) def func (w ):
return result .table [1].data [-1]+step *lf .fun (t +step ,w )-w
ynext =ynextinit ;h =step /10.
while func (ynext )>1.e-5:
ynext -=func (ynext )/((func (ynext +h )-func (ynext -h ))/2./h )
result .append ([t + step , ynext ])
t += step
return result
预报矫正法股份转让合同
预报矫正法也叫显⽰梯形法，预报矫正法的思想是先由当前点的值求出下⼀点的预报值：
再⽤当前点的值和预报点的值矫正下⼀个点的预报值：
也可以写作：
{w =y 00w =w +hf t ,w i +1i (i +1i +1)
(3)w i +1w i w i +1′w =i +1′w +i hf x ,w (i i )(3)
{w =w +h i +1i 2f t ,w +f t +h ,w (i i )(i i +1′)w =y 00(3)
{w =y 00w =w +f t ,w +f t +h ,w +hf t ,w i +1i 2h ((i i )(i i (i i )))(4)
def forecastCorrection (lf , tend , step =None ):
if step is None :
step = (tend - lf .t0) / 1e3
from ArrayTable import ArrayTable
result = ArrayTable (2, 0)
result .tTableHeader (["$t$", "y(t) solved with forecast correction"])
result .tTableUnit (["/", "/"])
苹果悬浮球result .append ([lf .t0, lf .y0])
t = lf .t0
while t < tend :
dydt =lf .fun (t , result .table [1].data [-1])
ynextpridict =result .table [1].data [-1] + dydt * step
ynext =result .table [1].data [-1]+step /2.*(dydt +lf .fun (t +step ,ynextpridict ))
result .append ([t + step , ynext ])
t += step
return result
4阶龙格-库塔⽅法(RK4)
def Rungekutta4(lf , tend , step =None ):
if step is None :
step = (tend - lf .t0) / 1e3
from ArrayTable import ArrayTable
result = ArrayTable (2, 0)
result .tTableHeader (["$t$", "y(t) solved with runge kutta 4"])
不可以的英文
result .tTableUnit (["/", "/"])
result .append ([lf .t0, lf .y0])
t = lf .t0
while t < tend :
s1=lf .fun (t ,result .table [1].data [-1])
s2=lf .fun (t +step /2.,result .table [1].data [-1]+step /2.*s1)
s3=lf .fun (t +step /2.,result .table [1].data [-1]+step /2.*s2)
s4=lf .fun (t +step ,result .table [1].data [-1]+step *s3)
ynext =result .table [1].data [-1]+step /6.*(s1+2*s2+2*s3+s4)
压强的物理意义result .append ([t + step , ynext ])
t += step
return result
结果展⽰
下图展⽰了使⽤不同⽅法求解
的结果，其中隐式欧拉⽅采⽤的步长为0.01，由于步长较⼤的原因，导致全局截断误差较⼤。常微分⽅程与最优化问题的转化
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧w =w +s +2s +2s +s i +1i 6h
(1234)s =f t ,w 1(i i )s =f t +,w +s 2(i 2h i 2h 1)s =f t +,w +s 3(i 2h i 2h 2)
s =f t +h ,w +hs 4(i i 3)
w =y 00{=ye dt
dy t
y 0=1()
变分法
先将⽅程齐次化，令，则初值问题化为：
化为变分问题：
令求泛函问题：找⼀个试探函数空间，其基函数为神经⽹络解法1. z (t )=y (t )−y 0{=f t ,z +y ,t ∈a ,b dt
dz (0)[]z a =0()−f t ,z +y vdt ∫a b (dt dz (0))=zv ∣−z dt −f t ,z +y vdt a b ∫a b dt dv ∫a b (0)=z b v b −z dt −f t ,z +y vdt =0
()()∫a b dt dv ∫a b (0)J (u )=u b −()2u dt −∫a b dt du f t ,u +y udt
∫a b (0){u ∈S ,S =v ∣v ∈C a ,b ,v a =00101
{1[]()}J u ≤J w ,∀w ∈S ()()01φ(x )=i (a −t )t ,i =i 0,1,2,...
