Matlab微分方程求解

引言

 数模课上讲的微分方程部分的知识，在之前一直是笔者弱项，在此将基础知识稍作整理后汇总成此博文，不足之处，望笔者多多指正。

微分方程求解

Matlab求解析解

 在Matlab中解析解常用函数库dslove来进行求解，常用的语法格式：

dsolve(‘方程1’,‘方程2’,…,‘方程n’,‘初始条件’,‘自变量’)
说明：在表达微分方程时，用字母D表示求微分，D2、D3等
表示求高阶微分.任何D后所跟的字母为因变量，自变量可以指
定或由系统规则选定为缺省

• 例子1：求取$\frac{du}{dt}=1+u^2$的通解。

dsolve('Du=1+u^2','t')

• 例子2：求取$$\{\begin{array}{l}\frac{d^{2}y}{d{x}^2}+4\frac{dy}{dx}+29y=0\\ y(0)=0,y^{\prime }(0)=15\end{array}$$ 的通解

y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')

• 例子3：求方程 $$\{\begin{array}{l}\frac{dx}{dt}=2x-3y+3z\\ \frac{dy}{dt}=4x-5y+3z\\ \frac{dz}{dt}=4x-4y+2z\end{array}$$通解

x,y,z]=dsolve('Dx=2*x-3*y+3*z',...
　　'Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't')

Matlab求取数值解

 在生产和科研中所处理的微分方程往往很复杂,且大多得不出一般解．而实际中的对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式。

求解原理

• 欧拉法
   在求解数值解过程中假设$x_{i+1}-x_i=h,i=0,1,2,\cdots ,n-1$,那么可以尝试用离散的方式去接微分方程：$$\{\begin{array}{l}{y}^{\prime }=f(x,y)\\ y\left(x_0\right)=y_0\end{array}$$
  如果h步长较小，有：$$y^{\prime }(x)\approx \frac{y(x+h)-y(x)}{h}$$
  可以得等到对应的欧拉公式：$$\{\begin{array}{l}{y}_{i+1}=y_i+hf\left(x_i,y_i\right)\\ y_0=y\left(x_0\right)\end{array}i=0,1,2,\cdots ,n-1$$
• 使用数值积分(改进欧拉)

对方程f(x,y)两边由$x_i$到$x_{i+1}$进行积分，并利用梯形公式有：$$\begin{array}{rlr}y\left(x_{i+1}\right)-y\left(x_i\right)& ={\int }_{x_i}^{x_{i+1}}f(t,y(t))dt& \approx \frac{x_{i+1}-x_i}{2}\left[f\left(x_i,y\left(x_i\right)\right)+f\left(x_{i+1},y\left(x_{i+1}\right)\right)\right]\end{array}$$
得到：$$\{\begin{array}{l}{y}_{i+1}=y_i+\frac{h}{2}\left[f\left(x_i,y_i\right)+f\left(x_{i+1},y_{i+1}\right)\right]\\ y_0=y\left(x_0\right)\end{array}$$
结合欧拉公式使用可以得到：
$$\{\begin{array}{l}{y}_{i+1}^{(0)}=y_i+hf\left(x_i,y_i\right)\\ y_{i+1}^{(k+1)}=y_i+\frac{h}{2}\left[f\left(x_i,y_i\right)+f\left(x_{i+1},y_{i+1}^{(k)}\right)\right]k=0,1,2\end{array}$$
满足精度条件后继续下一步的计算。

• 此外还可以结合泰勒公式继续得到比如龙格—库塔法、线性多步法等。

matlab求解

 语法格式：

[t，x]=solver(’f’,ts,x0,options)
(偷懒直接截图)

• 例子4：$$\{\begin{array}{c}\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}-1000\left(1-x^2\right)\frac{\mathrm{d}x}{\mathrm{d}t}-x=0\\ x(0)=2;x^{\prime }(0)=0\end{array}$$
  求解程序：

%例子4函数
function dy=vdp1000(t,y)
    dy=zeros(2,1);
    dy(1)=y(2);
    dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
end

命令窗口：
%例子4
t0=0;
tf=3000;
[T,Y]=ode15s('vdp1000',[0 3000],[2 0]); 
     plot(T,Y(:,1),'-')

运行结果：

• 例子5：$$\{\begin{array}{c}{y}_1^{\prime }=y_2{y}_3\\ y_2^{\prime }=-y_1{y}_3\\ y_3^{\prime }=-0.51y_1{y}_2\\ y_1(0)=0,y_2(0)=1,y_3(0)=1\end{array}$$
  求解程序：

%求解例子5
function dy=rigid(t,y)
       dy=zeros(3,1);
       dy(1)=y(2)*y(3);
       dy(2)=-y(1)*y(3);
       dy(3)=-0.51*y(1)*y(2);
end

命令行：
t0=0;
 tf=12;
 [T,Y]=ode45('rigid',[0 12],[0 1 1]);
     plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')

结果：

小结

 在笔者本科阶段的学习过程中，微分求解方程的学习经历课程主要有：常微分方程、数值分析、数学建模。希望在未来的学习过程中将其进一步的串成网络。最后本文不足之处望多多指正。

一段

某程序员夫妇新婚，一年之后喜得贵子，取名"灵灵"
过一年后又喜得一女，取名"灵伊"
两年之后得子"伊灵"
两年之后，夫妇商定为得圆满再生一子，取名"伊伊"
不料产科发现所怀为双胞胎，夫欲减胎，妻不允，冥思许久，对夫曰:“老五就叫’忆初’吧…”