【Matlab】一、解常微分方程ODE
文章目录
- 求解常微分方程 ODE
- （1）求解解析解
- （2）求解数值解
求解常微分方程 ODE
在matlab中，我们可以求解常微分方程的解析解，和数值解，一般使用dsolve来求解常微分方程的解析解，使用类似于ode45的求解器来求解常微分方程的数值解。
（1）求解解析解
求解解析解，例如求解该方程的解析解
dydx=3x2+1\frac{dy}{dx}=3x^2 + 1 dxdy=3x2+1
只需要在命令行中输入
dsolve('Dy=3*x^2+1', 'x')
或者是加上初始条件，求该方程在该初始条件下的解
dydx=3x2,y∣x=0=2\frac{dy}{dx}=3x^2, y|_{x=0}=2 dxdy=3x2,y∣x=0=2
在命令行输入
dsolve("Dy = 3*x^2", "y(0)=2", "x")
例如求一个常微分方程组
{x˙=y,y¨−y˙=0,x∣t=0=1;y˙∣t=0=1\left\{ \begin{array}{lr} \dot{x}=y, \\ \ddot{y}-\dot{y}=0, \quad x|_{t=0}=1; \dot{y}|_{t=0}=1 \end{array} \right. {x˙=y,y¨−y˙=0,x∣t=0=1;y˙∣t=0=1
在命令行输入
[x, y] = dsolve('Dx=y, D2y-Dy=0', 'x(0)=2, y(0)=2, Dy(0)=1')
（2）求解数值解
求数值解，有一些非线性的常微分方程是不能求出解析解的，我们一般求取其在一段区间内的数值解，采用迭代的方式来求解数值解，ode是matlab专门用于解微分方程的功能函数，具体的说明如下：
ode45是解决问题的首选，如果长时间没有结果，那么则采用ode15s试试。下面介绍ode45的函数格式
%函数格式
%[T,Y] = ode45(‘odefun’,tspan,y0)
%[T,Y] = ode45(‘odefun’,tspan,y0,options)
%[T,Y,TE,YE,IE] = ode45(‘odefun’,tspan,y0,options)
%sol = ode45(‘odefun’,[t0 tf],y0...)
%其中: odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名；
% tspan 是求解区间 [t0 tf]，或者一系列散点[t0,t1,...,tf]；
% y0 是初始值向量
% T 返回列向量的时间点
% Y 返回对应T的求解列向量
% options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等
% TE 事件发生时间
% YE 事件发生时之答案
% IE 事件函数消失时之指针i
首先得将微分方程标准化，类似于选择状态变量写状态空间表达式，对于一般的微分方程
{F(y,y˙,y¨,…,y(n−1),t)=0y(t0)=c0,y˙(t)=c1,…,y(n−1)(t0)=cn−1\left\{ \begin{array}{lr} F(y,\dot{y},\ddot{y},\dots,y^{(n-1)},t)=0 \\ y(t_0)=c_0, \dot{y}(t)=c_1, \dots, y^{(n-1)}(t_0)=c_{n-1} \end{array} \right. {F(y,y˙,y¨,…,y(n−1),t)=0y(t0)=c0,y˙(t)=c1,…,y(n−1)(t0)=cn−1
首先选择选择一组微分作为状态变量
x=[x1x2x3⋮xn]=[yy˙y¨⋮y(n−1)]\bold x = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \\x_{n} \end{bmatrix} = \begin{bmatrix} y \\ \dot{y} \\ \ddot{y} \\ \vdots \\y^{(n-1)} \end{bmatrix} x=⎣⎢⎢⎢⎢⎢⎡x1x2x3⋮xn⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡yy˙y¨⋮y(n−1)⎦⎥⎥⎥⎥⎥⎤
然后将待求解的微分方程F(y,y˙,y¨,…,y(n−1),yn,t)=0F(y,\dot{y},\ddot{y},\dots,y^{(n-1)},y^n,t)=0F(y,y˙,y¨,…,y(n−1),yn,t)=0，写成yn=G(y,y˙,y¨…,yn−1,t)=G(x1,x2,x3…,xn,t)y^{n}=G(y,\dot{y},\ddot{y} \dots,y^{n-1}, t)=G(x_1,x_2,x_3 \dots,x_{n},t)yn=G(y,y˙,y¨…,yn−1,t)=G(x1,x2,x3…,xn,t)的形式，然后写出x˙\dot{\bold x}x˙
x˙=[x1˙x2˙x3˙⋮xn˙]=[y˙y¨⋮y(n−1)y(n)]=[x2x3x4⋮G(x1,x2,x3…,xn,t)]\bold{\dot{x}} = \begin{bmatrix} \dot{x_1} \\ \dot{x_2} \\ \dot{x_3} \\ \vdots \\ \dot{x_n} \end{bmatrix} = \begin{bmatrix} \dot{y} \\ \ddot{y} \\ \vdots \\ y^{(n-1)} \\ y^{(n)} \end{bmatrix} = \begin{bmatrix} x_2 \\ x_3 \\ x_4 \\ \vdots \\ G(x_1,x_2,x_3 \dots,x_{n},t) \end{bmatrix} x˙=⎣⎢⎢⎢⎢⎢⎡x1˙x2˙x3˙⋮xn˙⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡y˙y¨⋮y(n−1)y(n)⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡x2x3x4⋮G(x1,x2,x3…,xn,t)⎦⎥⎥⎥⎥⎥⎤
上述步骤便是我们需要在odefun中完成的，举一个示例：在时间区间t=[3.9,4]t=[3.9,4]t=[3.9,4]，求解微分方程
y′′=−ty+ety′+3sin2t,y∣t=3=8,y′∣t=3=2y''=-ty+e^ty'+3sin2t, y|_{t=3}=8, y'|_{t=3}=2 y′′=−ty+ety′+3sin2t,y∣t=3=8,y′∣t=3=2
那么即
x˙=[x[1]−t∗x[1]+et∗x[2]+3sin2t]\bold{\dot{x}} = \begin{bmatrix} \bold{x}[1] \\ -t*\bold{x}[1]+e^t*\bold{x}[2]+3sin2t \end{bmatrix} x˙=[x[1]−t∗x[1]+et∗x[2]+3sin2t]
%在odefun.m脚本文件中完成以下内容
function dxdt = odefun(t, x)dxdt = zeros(2, 1); %初始化为 2 x 1的零矩阵dxdt(1) = x(2);dxdt(2) = -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
%在main.m脚本文件中完成以下内容
tspan = [3.9, 4];
y0 = [8, 2];
[t, y] = ode45('odefun', tspan, y0); %x的第一列为y，第二列为y’。如果遇到变量不是列向量形式的，可以考虑利用reshape函数做矩阵变换。
plot(t, y(:,1), '-o', t, y(:,2), '-*');
legend('y', "y'");
xlabel('t');
