Matlab求常微分方程组的数值解

上篇博客介绍了Matlab求解常微分方程组解析解的方法：博客地址
微分方程组复杂时，无法求出解析解时，就需要求其数值解，这里来介绍。
以下内容按照Matlab官方文档提供的方程来展开（提议多看官方文档）

介绍一下核心函数ode45()

一般形式：[t,y] = ode45(odefun,tspan,y0) 其中 tspan = [t0 tf]
功能介绍：求微分方程组 y′=f(t,y) 从 t0 到 tf 的积分，初始条件为 y0。解数组 y 中的每一行都与列向量 t 中返回的值相对应。

1. 一阶微分方程求解(简单调用即可)

方程：y’=2*t
代码：

tspan=[1 6]; %定义自变量x的取值空间为1-6
y0=0;%定义因变量的初值，当x=1（x取值空间的第一个数）时，y0=0
[t,y]=ode45(@(t,y) 2*t,tspan,y0); %定义函数y'=2*t,使用ode45求解
plot(t,y,'-o'); %绘制求得的数值曲线

说明：简单的odefun参数就是这个形式，@(x,y) fun 中间是空格，不能用逗号
结果：

2. 二阶微分方程求解(引入函数文件)

方程：范德波尔方程 y1’’-u(1-y1²)*y1’+y1=0;（这里设u=1）
代码：
定义输入的方程，以函数形式定义

function dydt=odefun(t,y)
%二阶方程为y1''-(1-y1²)*y1'+y1=0;
%降阶为两个方程：y1'=y2; 
%               y2'=(1-y1²)*y2-y1;
%t虽然没有使用，但必须要作为参数写入
dydt=[y(2);(1-y(1)^2)*y(2)-y(1)];

求解作图

tspan=[0 20]; %定义自变量x的取值空间为0-20
y0=[2;0];%定义因变量的初值，当x=0时，y1=2,y2=y1'=0;
[t,y]=ode45(@odefun,tspan,y0); %使用ode45求解
%%下面为作图过程，不解释
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

说明：先进行降阶，使y2=y1’，这时候y2’=y1’’，然后把两个方程写成列向量形式就行了
结果：

3. 求解微分方程组（和2类似）

这里就和求解二阶方程类似的，只不过不需要降阶，仍旧需要一个函数来定义方程组。我们这里不用官方文档的例子，用同学的循坏摆问题来进行演示。
方程：
给定的初值（w接近0，但实际上不能设置为0）：

代码：
定义输入的方程

function dRvw=func(t,Rvw)
%% 函数功能：为ode45提供微分方程
%输入：t:时间序列，就是θ；Rvw:因变量,Rvw(1)代表R,Rvw(2)代表v,Rvw(3)代表w
%输出：dRvw:因变量的一阶微分,dRvw(1)代表dR,dRvw(2)代表dv,dRvw(3)代表dw
%% 初始化因变量的一阶微分，3×1的向量
dRvw=zeros(3,1);
%% 参数初始化
r=0.01;u=0.1;g=9.8;M=10;m=1;
%% 输入微分方程式
dRvw(1)=-Rvw(2)/Rvw(3)-r;
dRvw(2)=(M*g-(m*Rvw(3)^2*Rvw(1)+m*g*sin(t)) *exp(u*(t+pi/2)))/(Rvw(3)*(M+m));
dRvw(3)=Rvw(3)*r/Rvw(1)+2*Rvw(2)/Rvw(1)+g*cos(t)/(Rvw(3)*Rvw(1));

求解作图

%% 参数初始化 start_Theta是θ的初始值 end_Theta是θ的结束值
%R是半径初值；v是线速度初值；w是角速度初值
start_Theta=0;end_Theta=2*pi;R=1;v=0;w=1e-5;
%% 使用ode45方法计算微分方程组func的数值解
%func是带有方程组的函数 
%[start_Theta end_Theta]是自变量范围
%[R;v;w]是方程初值
%T是自变量的数组，Rvw是对应的因变量的数值。Rvw(:,1)=R;Rvw(:,2)=v;Rvw(:,3)=w;
[T,Rvw]=ode45(@func,[start_Theta end_Theta],[R;v;w]);
%% 组合自变量和因变量。TRvw(:,1)=θ;TRvw(:,2)=R;TRvw(:,3)=v;Trvw(:,4);
TRvw=[T,Rvw];
%% 绘制自变量-因变量图像.figure1是R-θ,figure2是v-θ,figure3是w-θ
plot(T,Rvw(:,1),'-o',T,Rvw(:,2),'-o',T,Rvw(:,3),'-o')
title('自变量-因变量图像');
xlabel('Angle θ');
ylabel('Solution');
legend('R','v','w')

说明：注释的应该是比较清楚的，把三个方程写成列向量的形式就行

PS:有些人和我说不能运行，然后我看了他们出错的原因，有点儿哭笑不得。出错的基本上都是运行上面的dRvw=func（t，Rvw）这个函数的。说明一下，这是有参数的函数，不给参数不能直接运行的。下面的求解作图脚本才是需要运行的哈，它调用了函数，才得到的结果。如果大家还发现什么问题，欢迎私戳或评论。

PS+ 有了PS之后，还是很多人问我参数的问题，我在这里直接把文件给大家：cupt.zip
提取码：6k8n。
解压后，有两个文件，大家直接运行cupt.m即可。

结果：

4. 更多形式

讲到这里，大部分我们用到的微分方程形式都可以求解了，Matlab还支持带有时变项和额外参数的微分方程求解，这里不再赘述，大家可以自行参阅官方文档。