IVP问题

对于常微分方程的初值问题
{dydt=f(t,y),t∈[a,b]y(a)=y0(1)\left\{ \begin{array}{l} \frac{{dy}}{{dt}} = f(t,y),t \in \left[ {a,b} \right]\\ y\left( a \right) = {y_0} \end{array} \right.\tag{1} {dtdy=f(t,y),t∈[a,b]y(a)=y0(1)
有以下解的存在性、唯一性定理：

假设f(t,y)f(t,y)f(t,y)定义在集合[a,b]×[α,β][a,b]\times[\alpha,\beta][a,b]×[α,β]并且初值α<y0<β\alpha<y_0<\betaα<y0<β，函数对于变量yyy是Lipschitz连续，则在aaa与bbb之间存在ccc，使得初值问题
{dydt=f(t,y),t∈[a,b]y(a)=y0t∈[a,c](2)\left\{ \begin{array}{l} \frac{{dy}}{{dt}} = f(t,y),t \in \left[ {a,b} \right]\\ y\left( a \right) = {y_0}\\ t\in[a,c] \end{array} \right.\tag{2} ⎩⎨⎧dtdy=f(t,y),t∈[a,b]y(a)=y0t∈[a,c](2)
有唯一解y(t)y(t)y(t)。而且，如果fff在[a,b]×[−∞,+∞][a,b]\times[-\infty,+\infty][a,b]×[−∞,+∞]上是Lipschitz连续，则其在[a,b][a,b][a,b]上存在唯一解。

定义一个ODEsolver类来使用各种方法求解常微分方程，使用待求解的函数和初始值来初始化该类：

class ODEsolver:def __init__(self,fun,t0,y0):self.fun=funself.t0=t0self.y0=y0

Euler显式格式

Euler方法为：
{w0=y0wi+1=wi+hf(ti,wi)(2)\left\{ \begin{array}{l} {w_0} = {y_0}\\ {w_{i + 1}} = {w_i} + hf\left( {{t_i},{w_i}} \right) \end{array} \right.\tag{2} {w0=y0wi+1=wi+hf(ti,wi)(2)
代码为：

    def eulersolver(self,tend,step=None):if step is None:step=(tend-self.t0)/1e3from ArrayTable import ArrayTableresult=ArrayTable(2,0)result.setTableHeader(["$t$","$y(t)$"])result.setTableUnit(["/","/"])result.append([self.t0,self.y0])t=self.t0while t<tend:ynext=result.table[1].data[-1]+self.fun(t,result.table[1].data[-1])*stepresult.append([t+step,ynext])t+=stepreturn result

隐式Euler法

隐式Euler法的原理为：
{w0=y0wi+1=wi+hf(ti+1,wi+1)(3)\left\{ \begin{array}{l} {w_0} = {y_0}\\ {w_{i + 1}} = {w_i} + hf\left( {{t_{i + 1}},{w_{i + 1}}} \right) \end{array} \right.\tag{3} {w0=y0wi+1=wi+hf(ti+1,wi+1)(3)
迭代开始前首先要预估一个wi+1w_{i+1}wi+1的初始值，一般有两种选择，使用wiw_iwi作为初始值或者用显示方法先预估一个值wi+1′w_{i+1}'wi+1′，这里用预报矫正法来估计其值。

    def implicitEuler(self,tend, step=None):if step is None:step = (tend - self.t0) / 1e3from ArrayTable import ArrayTableresult = ArrayTable(2, 0)result.setTableHeader(["$t$", "$y(t)$"])result.setTableUnit(["/", "/"])result.append([self.t0, self.y0])t = self.t0while t<tend:temp=self.fun(t,result.table[1].data[-1])ynextinit=result.table[1].data[-1]+step/2.*(temp+self.fun(t+step,result.table[1].data[-1] + temp * step))def func(w):return result.table[1].data[-1]+step*self.fun(t+step,w)-wynext=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 += stepreturn result

预报矫正法

预报矫正法也叫显示梯形法¹，预报矫正法的思想是先由当前点的值求出下一点的预报值：
wi+1′=wi+hf(xi,wi)(3)w_{i + 1}' = {w_i} + hf\left( {{x_i},{w_i}} \right)\tag{3} wi+1′=wi+hf(xi,wi)(3)
再用当前点的值和预报点的值矫正下一个点的预报值：
{wi+1=wi+hf(ti,wi)+f(ti+h,wi+1′)2w0=y0(3)\left\{ \begin{array}{l} {w_{i + 1}} = {w_i} + h\frac{{f\left( {{t_i},{w_i}} \right) + f\left( {{t_i} + h,w_{i + 1}'} \right)}}{2}\\ {w_0} = {y_0} \end{array} \right.\tag{3} {wi+1=wi+h2f(ti,wi)+f(ti+h,wi+1′)w0=y0(3)
也可以写作：
{w0=y0wi+1=wi+h2(f(ti,wi)+f(ti+h,wi+hf(ti,wi)))(4)\left\{ \begin{array}{l} {w_0} = {y_0}\\ {w_{i + 1}} = {w_i} + \frac{h}{2}\left( {f\left( {{t_i},{w_i}} \right) + f\left( {{t_i} + h,{w_i} + hf\left( {{t_i},{w_i}} \right)} \right)} \right) \end{array} \right.\tag{4} {w0=y0wi+1=wi+2h(f(ti,wi)+f(ti+h,wi+hf(ti,wi)))(4)

    def forecastCorrection(self, tend, step=None):if step is None:step = (tend - self.t0) / 1e3from ArrayTable import ArrayTableresult = ArrayTable(2, 0)result.setTableHeader(["$t$", "y(t) solved with forecast correction"])result.setTableUnit(["/", "/"])result.append([self.t0, self.y0])t = self.t0while t < tend:dydt=self.fun(t, result.table[1].data[-1])ynextpridict=result.table[1].data[-1] + dydt * stepynext=result.table[1].data[-1]+step/2.*(dydt+self.fun(t+step,ynextpridict))result.append([t + step, ynext])t += stepreturn result

4阶龙格-库塔方法(RK4)

{wi+1=wi+h6(s1+2s2+2s3+s4)s1=f(ti,wi)s2=f(ti+h2,wi+h2s1)s3=f(ti+h2,wi+h2s2)s4=f(ti+h,wi+hs3)w0=y0\left\{ \begin{array}{l} {w_{i + 1}} = {w_i} + \frac{h}{6}\left( {{s_1} + 2{s_2} + 2{s_3} + {s_4}} \right)\\ {s_1} = f\left( {{t_i},{w_i}} \right)\\ {s_2} = f\left( {{t_i} + \frac{h}{2},{w_i} + \frac{h}{2}{s_1}} \right)\\ {s_3} = f\left( {{t_i} + \frac{h}{2},{w_i} + \frac{h}{2}{s_2}} \right)\\ {s_4} = f\left( {{t_i} + h,{w_i} + h{s_3}} \right)\\ {w_0} = {y_0} \end{array} \right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧wi+1=wi+6h(s1+2s2+2s3+s4)s1=f(ti,wi)s2=f(ti+2h,wi+2hs1)s3=f(ti+2h,wi+2hs2)s4=f(ti+h,wi+hs3)w0=y0

    def Rungekutta4(self, tend, step=None):if step is None:step = (tend - self.t0) / 1e3from ArrayTable import ArrayTableresult = ArrayTable(2, 0)result.setTableHeader(["$t$", "y(t) solved with runge kutta 4"])result.setTableUnit(["/", "/"])result.append([self.t0, self.y0])t = self.t0while t < tend:s1=self.fun(t,result.table[1].data[-1])s2=self.fun(t+step/2.,result.table[1].data[-1]+step/2.*s1)s3=self.fun(t+step/2.,result.table[1].data[-1]+step/2.*s2)s4=self.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 += stepreturn result

结果展示

下图展示了使用不同方法求解
{dydt=yety(0)=1\left\{ \begin{array}{l} \frac{{dy}}{{dt}} = y{e^t}\\ y\left( 0 \right) = 1 \end{array} \right. {dtdy=yety(0)=1
的结果，其中隐式欧拉方采用的步长为0.01，由于步长较大的原因，导致全局截断误差较大。

常微分方程与最优化问题的转化

变分法

先将方程齐次化，令z(t)=y(t)−y0z(t)=y(t)-y_0z(t)=y(t)−y0，则初值问题化为：
{dzdt=f(t,z+y0),t∈[a,b]z(a)=0\left\{ \begin{array}{l} \frac{{dz}}{{dt}} = f\left( {t,z + {y_0}} \right),t \in \left[ {a,b} \right]\\ z\left( a \right) = 0 \end{array} \right. {dtdz=f(t,z+y0),t∈[a,b]z(a)=0
化为变分问题：
∫ab(dzdt−f(t,z+y0))vdt=zv∣ab−∫abzdvdtdt−∫abf(t,z+y0)vdt=z(b)v(b)−∫abzdvdtdt−∫abf(t,z+y0)vdt=0\begin{array}{l} \int_a^b {\left( {\frac{{dz}}{{dt}} - f\left( {t,z + {y_0}} \right)} \right)vdt} \\ = zv|_a^b - \int_a^b {z\frac{{dv}}{{dt}}dt - \int_a^b {f\left( {t,z + {y_0}} \right)vdt} } \\ = z\left( b \right)v\left( b \right) - \int_a^b {z\frac{{dv}}{{dt}}dt} - \int_a^b {f\left( {t,z + {y_0}} \right)vdt} = 0 \end{array} ∫ab(dtdz−f(t,z+y0))vdt=zv∣ab−∫abzdtdvdt−∫abf(t,z+y0)vdt=z(b)v(b)−∫abzdtdvdt−∫abf(t,z+y0)vdt=0
令J(u)=u(b)2−∫abududtdt−∫abf(t,u+y0)udtJ(u)=u\left( b \right)^2 - \int_a^b {u\frac{{du}}{{dt}}dt} - \int_a^b {f\left( {t,u + {y_0}} \right)udt}J(u)=u(b)2−∫abudtdudt−∫abf(t,u+y0)udt
求泛函问题：
{u∈S01,S01={v∣v∈C1[a,b],v(a)=0}J(u)≤J(w),∀w∈S01\left\{ \begin{array}{l} u \in S_0^1,S_0^1 = \left\{ {v|v \in {C^1}\left[ {a,b} \right],v\left( a \right) = 0} \right\}\\ J\left( u \right) \le J\left( w \right),\forall w \in S_0^1 \end{array} \right. {u∈S01,S01={v∣v∈C1[a,b],v(a)=0}J(u)≤J(w),∀w∈S01
找一个试探函数空间，其基函数为φi(x)=(a−t)ti,i=0,1,2,...\varphi_i(x)=(a-t)t^i,i=0,1,2,...φi(x)=(a−t)ti,i=0,1,2,...

神经网络解法

Timothy Sauer.数值分析第2版中文版[M].北京：机械工业出版社.2020. ↩︎

常微分方程求解器ODE solver相关推荐

【工具笔记】Microsoft数学求解器Math Solver
[工具笔记]Microsoft数学求解器Math Solver 工具笔记用于记录各种有用的工具,这里记录的是一个由Microsoft提供的数学求解器Math Solver. 可以用于求解代数,三角学, ...
matlab中solver函数_Simulink求解器（Solver）相关知识
更多精彩内容参见专业MATLAB技术交流平台--MATLAB技术论坛http://www.matlabsky.com 1.变步长(Variable-Step)求解器可以选择的变步长求解器有:ode4 ...
JAX-FLUIDS：可压缩两相流的完全可微高阶计算流体动力学求解器
原文来自微信公众号"编程语言Lab":论文精读 | JAX-FLUIDS:可压缩两相流的完全可微高阶计算流体动力学求解器搜索关注"编程语言Lab"公众号(HW ...
CST入门——求解器简介与时域、频域和积分求解器设置
目录 1. 高频电磁仿真求解器简介 1.1. 时域求解器 Time Domain Solver(主) 1.2. 频域求解器 Frequency Domain Solver(主) 1.3. 本征模求解器 ...
【JY】No.7.1力学架构结构力学求解器(SM)使用教程
软件讲解示例之理论(电算/手算)分析思路: 1.结构建模,并完成结构假定,如定义梁端弯矩释放(详见第二章). 2.线性屈曲分析(是否失稳破坏分析): 3.强度/挠度计算:(构件本身强度是否达标) 4. ...
CST微波工作室学习笔记—9.求解器
CST微波工作室_求解器详解 - 求解器的分类 Time Domain Solver--时域求解器 Frequency Domain Solver--频域求解器 Eigenmode Solver--本 ...
【AI】求解器SGD、BGD、MBGD等详解
参考博客: ***** 深度学习必备:随机梯度下降(SGD)优化算法及可视化: **** 深度学习--优化器算法Optimizer详解(BGD.SGD.MBGD.Momentum.NAG.Adagra ...
matlab中离散数值求解器在哪,matlab - 在ode MATLAB求解器上使用中间值 - SO中文参考 - www.soinside.com...
我正在使用刚性求解器(ode15s)对ODE系统进行时间积分.它工作正常,但我想加快速度. 方程组以状态空间形式给出: function [dx] = fun(t,x,M,C,K,other_para ...
MATLAB 数学应用微分方程常微分方程求解非刚性ODE
本文介绍两个使用 ode45 来求解非刚性常微分方程的示例.MATLAB拥有三个非刚性 ODE 求解器. ode45 ode23 ode113 对于大多数非刚性问题,ode45 的性能最佳.但对于允许 ...
关于求解微分方程——初学Matlab里的 ODE求解器
学习背景最近想挖掘一下自己项目的理论深度,于是找到了老师.在老师的建议下,我们开始了漫长的研读老师的论文的旅程(论文名:Optimal Design of Adaptive Robust Contr ...

常微分方程求解器ODE solver

文章目录