微分方程的解法（Matlab）

微分方程分为刚性微分方程和非刚性微分方程，在数值解法中的表现和行为特性上存在显著差异。

刚性微分方程（Stiffness Equation）是指其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。这类方程通常包含多个时间尺度，且至少有一个时间尺度非常小，导致解的变化速率差异很大。
刚性微分方程对数值解法的微小误差非常敏感，传统的显式数值方法（如显式欧拉法）在求解时可能遇到稳定性问题，因为解的变化速率差异大，显式方法可能由于稳定性限制而需要非常小的步长，从而导致计算效率低下。
求解方法：需要使用隐式方法（如隐式欧拉法、后退欧拉法、亚当斯-穆森方法、BDF、SDIRK等）来求解，这些方法对时间步长的大小不那么敏感，能够在保证稳定性的同时允许使用较大的时间步长，从而提高计算效率。

非刚性微分方程

非刚性微分方程是指那些解的变化速率相对均匀，对数值解法误差不太敏感的系统。这类方程的所有时间尺度都相似，或者时间尺度的差异不大，解在整个时间范围内变化速度相对均匀。
非刚性微分方程可以使用较大的时间步长，并且显式数值方法通常足够稳定和有效。
求解方法：可以使用更广泛的数值方法进行求解，包括显式方法和隐式方法。显式方法（如显式欧拉法）因其简单性和易于实现，在求解非刚性微分方程时常常被采用。

Matlab数值解

非刚性常微分方程的解法

Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45，ode23， ode113，其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的首选方法，ode23 采用二三阶 RK 方法，ode113 采用的是多步法，效率一般比 ode45 高。

示例：

编写改进的 Euler 方法函数如下：

function [x,y]=eulerpro(fun,x0,xfinal,y0,n) 
if nargin<5,n=50;end 
h=(xfinal-x0)/n; 
x(1)=x0;y(1)=y0; 
for i=1:n 
    x(i+1)=x(i)+h; 
    y1=y(i)+h*feval(fun,x(i),y(i)); 
    y2=y(i)+h*feval(fun,x(i+1),y1); 
    y(i+1)=(y1+y2)/2; 
end

需要五个输入参数：fun 是微分方程 f(x,y) 的函数句柄，x0 是初始 x 值，xfinal 是求解区间的终点 x 值，y0 是初始 y 值，n 是将区间 [x0,xfinal] 分割成的小区间数（即步数），默认步长为50。函数返回两个数组 x 和 y，分别表示 x 值和对应的 y 值近似解。

function f=doty(x,y)
f=-2*y+2*x^2+2*x;

求解：

[x,y]=eulerpro('doty',0,0.5,1,10)

求得数值解如下：

ode23，ode45，ode113的使用

Matlab的函数形式是： [t，y]=solver('F'，tspan，y0)

这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程。

tspan=[t0，tfinal]是求解区间，y0是初值。

根据前边写好的函数文件：

[x,y]=ode45('doty',[0,0.5],1)

即可求出数值解：

或者用句柄函数也可以求出相同的数值解：

clc,clear;
f=@(x,y)-2*y+2*x^2+2*x;
[x,y]=ode45(f,[0,0.5],1)

刚性常微分方程的解法

Matlab的工具箱提供了几个解刚性常微分方程的功能函数，如ode15s，ode23s， ode23t，ode23tb，这些函数的使用同上述非刚性微分方程的功能函数。

高阶微分方程的解法

用 Matlab 解决此问题的函数形式为 [T,Y]=solver('F',tspan,y0) 这里 solver 为 ode45、ode23、ode113，输入参数 F 是用 M 文件定义的常微分方程组， tspan=[t0 tfinal]是求解区间，y0 是初值列向量。

把一阶方程组写成接受两个参数t 和 y ，返回一个列向量的 M 文件：

function dy=F(t,y); 
dy=[y(2);y(3);3*y(3)+y(2)*y(1)];

求上述常微分方程的数值解：

[T,Y]=ode45('F',[0 1],[0;1;-1])

这里 Y 和时刻 T 是一一对应的，Y(:,1)是初值问题的解，Y(:,2)是解的导数，Y(:,3)是解的二阶导数。

或者使用句柄函数可得相同结果：

F=@(t,y)[y(2);y(3);3*y(3)+y(2)*y(1)];
[T,Y]=ode45(F,[0 1],[0;1;-1])

示例：

书写 M 文件：

function dy=vdp1(t,y); 
dy=[y(2);(1-y(1)^2)*y(2)-y(1)];

调用 Matlab 函数。对于初值 y(0) = 2, y''(0) = 0 ，解为 :

[T,Y]=ode45('vdp1',[0 20],[2;0])

或使用句柄函数：

F=@(t,y)[y(2);(1-y(1)^2)*y(2)-y(1)];
[T,Y]=ode45(F,[0 20],[2;0])

将结果可视乎：

plot(T,Y(:,1),'-',T,Y(:,2),'--') 
title('Solution of van der Pol Equation,mu=1'); 
xlabel('time t'); 
ylabel('solution y'); 
legend('y1','y2');

van der Pol 方程， μ =1000 （刚性）

function dy=vdp1000(t,y); 
dy=[y(2);1000*(1-y(1)^2)*y(2)-y(1)];

得出结果，绘出结果如下：

[t,y]=ode15s('vdp1000',[0 3000],[2;0]); 
subplot(1,2,1)
plot(t,y(:,1),'o') 
title('Solution of van der Pol Equation,mu=1000'); 
xlabel('time t'); 
ylabel('solution y(:,1)');
subplot(1,2,2)
plot(t,y(:,2),'*')
title('Solution of van der Pol Equation,mu=1000'); 
xlabel('time t'); 
ylabel('solution y(:,2)');

常微分方程的解析解

在 Matlab 中，符号运算工具箱提供了功能强大的求解常微分方程的符号运算命令 dsolve。常微分方程在 Matlab 中按如下规定重新表达：

符号 D 表示对变量的求导。Dy 表示对变量 y 求一阶导数，当需要求变量的 n 阶导数时，用 Dn 表示，D4y 表示对变量 y 求 4 阶导数。

无初边值条件的常微分方程的解就是该方程的通解。

其使用格式为：

dsolve('diff_equation')

dsolve(' diff_equation'，'var')

式中 diff_equation 为待解的常微分方程，第 1 种格式将以变量 t 为自变量进行求解，第 2 种格式则需定义自变量 var。

syms x y 
diff_equ='x^2+y+(x-2*y)*Dy=0'; 
dsolve(diff_equ,'x')

求解带有初边值条件的常微分方程的使用格式为： dsolve('diff_equation'，'condition1，condition2，…'，'var') 其中 condition1，condition2，… 即为微分方程的初边值条件。

y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x')

求解常微分方程组的命令格式为：

dsolve('diff_equ1，diff_equ2，…'，'var')

dsolve('diff_equ1，diff_equ2，…'，'condition1，condition2，…'，'var')

第 1 种格式用于求解方程组的通解，第 2 种格式可以加上初边值条件，用于具体求解。

clc,clear 
equ1='D2f+3*g=sin(x)'; 
equ2='Dg+Df=cos(x)'; 
[general_f,general_g]=dsolve(equ1,equ2,'x') 
[f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

通解：

syms t 
a=[2,1,3;0,2,-1;0,0,2]; 
x0=[1;2;1]; 
x=expm(a*t)*x0

clc,clear 
syms t s 
a=[1,0,0;2,1,-2;3,2,1];ft=[0;0;exp(t)*cos(2*t)]; 
x0=[0;1;1]; 
x=expm(a*t)*x0+int(expm(a*(t-s))*subs(ft,s),s,0,t); 
x=simplify(x)