1、参考资料：

Solve system of differential equations - MATLAB dsolve- MathWorks 中国

Matlab求常微分方程组的解析解_y=520（2sinM-sin2M）的博客-CSDN博客_matlab求解常微分方程的解析解

Matlab求常微分方程组的解析解_y=520（2sinM-sin2M）的博客-CSDN博客_matlab求解常微分方程的解析解

Matlab求常微分方程组的数值解_y=520（2sinM-sin2M）的博客-CSDN博客_matlab求微分方程组的数值解

2、官方文档，YYDS

一阶：

二阶：

加上边界

求解微分方程组：

简 介：前面介绍了微分方程的解析解方法，同时也指出很多非线性微分方程是不存在解析解法的，需要使用数值解法对之进行研究。本文着重讨论基于 MATLAB/Simulink语言的各类微分方程的数值解方法。

关键词： 微分方程、数值解、MATLAB

§01 总述

一般微分方程的数值解法很大一类是关于微分方程初值问题的数值解法，这类问题需要用一阶显式的微分方程组来描述为

$$\dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t))$$

其中, ${\mathbf{x}}^{\mathrm{T}}(t)=\left[x_1(t),x_2(t),\cdots ,x_n(t)\right]$称为状态向量，${\mathbf{f}}^{\mathrm{T}}(\cdot )=\left[f_1(\cdot ),f_2(\cdot ),\cdots ,f_n(\cdot )\right]$可以是任意非线性函数。所谓初值问题是指，若已知初始状态${\mathbf{x}}_0={\left[x_1(0),\cdots ,x_n(0)\right]}^{\mathrm{T}}$，用数值求解方法求出在某个时间区间$t\in \left[0,t_{\mathrm{f}}\right]$内各个时刻状态变量 $\mathbf{x}(t)$ 的数值，这里 $t_{\mathrm{f}}$ 又称为终止时间。

其他类型微分方程求解必须先转换：

1. 单个高阶常微分方程处理方法

一个高阶常微分方程的一般形式为：

$$y^{(n)}=f\left(t,y,y^{\prime },\cdots ,y^{(n-1)}\right)$$

输出变量的各阶导数初始值为：

$$y(0),y^{\prime }(0),\cdots ,y^{(n-1)}(0)$$

选择一组状态变量：

$$x_1=y,x_2=y^{\prime },\cdots ,x_n=y^{(n-1)}$$

原高阶常微分方程模型变换为一阶微分方程组：

$$\{\begin{array}{rl}{x}_1^{\prime }& =x_2\\ x_2^{\prime }& =x_3\\ & ⋮\\ x_n^{\prime }& =f\left(t,x_1,x_2,\cdots ,x_n\right)\end{array}$$

初值为：

$$x_1(0)=y(0),x_2(0)=y^{\prime }(0),\cdots ,x_n(0)=y^{(n-1)}(0)$$

2. 高阶常微分方程组的变换方法

多元高阶常微分方程组的处理

$$\{\begin{array}{l}{x}^{(m)}=f\left(t,x,x^{\prime },\cdots ,x^{(m-1)},y,\cdots ,y^{(n-1)}\right)\\ \\ y^{(n)}=g\left(t,x,x^{\prime },\cdots ,x^{(m-1)},y,\cdots ,y^{(n-1)}\right)\end{array}$$

状态变量的选择不唯一，建议：选择如下状态变量

$$\begin{gathered} x_1=x,x_2=x^{\prime },\cdots ,x_m=x^{(m-1)}, \\ {x}_{m+1}=y,x_{m+2}=y^{\prime },\cdots ,x_{m+n}=y^{(n-1)} \end{gathered}$$

新的状态方程

$$\{\begin{array}{rl}{x}_1^{\prime }=& x_2\\ & ⋮\\ x_m^{\prime }=& f\left(t,x_1,x_2,\cdots ,x_{m+n}\right)\\ x_{m+1}^{\prime }& =x_{m+2}\\ & ⋮\\ x_{m+n}^{\prime }& =g\left(t,x_1,x_2,\cdots ,x_{m+n}\right)\end{array}$$

可以描述该方程，然后用 ode45 等求解。

§02 微分方程求解的误差与步长问题

1. 不能无限制地减小步长$h$的值的两条原因

• 减慢计算速度
• 增加累积误差

2. 在对微分方程求解过程中应采取的三个措施

• 选择适当的步长
• 改进近似算法精度
• 采用变步长方法

§03 函数调用格式（ode45）

[t, x] = ode45(fun,[t0, tf], x0)  % 直接求解
[t, x] = ode45(fun,[t0, tf], x0, options)  % 带有控制选项
[t, x] = ode45(fun,[t0, tf], x0, options, p1, p2, ...)  % 带有附加参数

当然，还存在其他求解函数，如ode15s、ode23等，不同的算法（函数）适合于不同类型的微分方程。

§04 微分方程的MATLAB描述

描述需要求解的微分方程组

  $\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}{\mathbf{x}}_{\mathrm{d}}=\mathrm{f}\mathrm{u}\mathrm{n}(\mathrm{t},\mathbf{x})$

  $\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}{\mathbf{x}}_{\mathrm{d}}=\mathrm{f}\mathrm{u}\mathrm{n}\left(\mathrm{t},\mathrm{x},{\mathrm{p}}_1,{\mathrm{p}}_2,\cdots \right)$

修改控制变量

options = odeset('RelTol', 1e-7) ;
options = odeset;  options.RelTol = 1e-7;

§05 应用举例：Lorenz方程

例1：无参数

题目： 求解下列Lorenz模型

$$\{\begin{array}{l}{\dot{x}}_1(t)=-\beta x_1(t)+x_2(t)x_3(t)\\ \\ {\dot{x}}_2(t)=-\rho x_2(t)+\rho x_3(t)\\ \\ {\dot{x}}_3(t)=-x_1(t)x_2(t)+\sigma x_2(t)-x_3(t)\end{array}$$

式中参数为 $\beta =8/3,\rho =10,\sigma =28$

初始条件为 $x_1(0)=x_2(0)=0,x_3(0)=ϵ,\mathrm{i}.\mathrm{e}.,ϵ=10^{-10}$

求解：

1. 写出标准型

$${\mathbf{x}}^{\prime }(t)=\left[\begin{array}{c}-\beta x_1(t)+x_2(t)x_3(t)\\ \\ -\rho x_2(t)+\rho x_3(t)\\ \\ -x_1(t)x_2(t)+\sigma x_2(t)-x_3(t)\end{array}\right],\mathbf{x}(0)=\left[\begin{array}{l}0\\ \\ 0\\ \\ ϵ\end{array}\right]$$

2. 微分方程的MATLAB描述

• M-文件描述

function y = lorenzeq(t, x)
	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];

• 匿名函数

f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];

3. 微分方程求解

• 匿名函数

t_final = 100; 
x0 = [0; 0; 1e-10]; 
f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
[t,x] = ode45(f, [0,t_final], x0); 
plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));

• M-文件描述

t_final = 100; 
x0 = [0; 0; 1e-10]; 
f = @(t,x)[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
[t,x] = ode45(@lorenzeq, [0,t_final], x0); 
plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3));

function y = lorenzeq(t, x)
	y = [-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3)];
end

例2：带有附加参数

引入附加参数的目的： 如果参数 $\beta$, $\rho$, $\sigma$ 改变，不需要修改原函数。

重新求解Lorenz方程

方式一

f = @(t,x,beta,rho,sigma)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
t_final=100; 
x0 = [0; 0; 1e-10]; 
b1 = 8/3; r1 = 10; s1 = 28; 
[t,x]=ode45(f, [0,t_final], x0, [], b1, r1, s1);

方式二

beta = 2; rho = 5; sigma = 20; 
f = @(t,x)[-beta*x(1) + x(2)*x(3); -rho*x(2) + rho*x(3); -x(1)*x(2) + sigma*x(2) - x(3)];
t_final = 100;
x0 = [0; 0; 1e-10]; 
[t2,x2] = ode45(f, [0,t_final], x0);

1、dsolve 函数

dsolve函数用于求常微分方程组的精确解，也称为常微分方程的符号解。如果没有初始条件或边界条件，则求出通解；如果有，则求出特解。

1)函数格式  

Y = dsolve(‘eq1,eq2,…’ , ’cond1,cond2,…’ , ’Name’)

其中，‘eq1,eq2,…’:表示微分方程或微分方程组;

            ’cond1,cond2,…’:表示初始条件或边界条件;

            ‘Name’:表示变量。没有指定变量时，matlab默认的变量为t；
 

利用dsolve函数进行求解。

 实例

下面呢，我要求解一个微分方程，现在求解一个物体的速度随时间变化的解析解。给一个初始速度为1的物体减速，减速的力与速度有关，还有一个恒力帮助减速。

微分方程：

初始条件：

%构建
syms y(t)
s= dsolve('-Dy=y +1','y(0) =1');
%作图
t = 0:0.1:5
y = eval(subs(s))
plot(t,y)

 在命令窗口行键入结果 's'，显示解析式：

        别林斯基曾经提到过，好的书籍是最贵重的珍宝。这句话看似简单，但其中的阴郁不禁让人深思. 要想清楚，微分方程，到底是一种怎么样的存在。 对我个人而言，微分方程不仅仅是一个重大的事件，还可能会改变我的人生。 现在，解决微分方程的问题，是非常非常重要的。 所以， 所谓微分方程，关键是微分方程需要如何写。 经过上述讨论， 微分方程的发生，到底需要如何做到，不微分方程的发生，又会如何产生。 微分方程，发生了会如何，不发生又会如何。 微分方程的发生，到底需要如何做到，不微分方程的发生，又会如何产生。 微分方程，发生了会如何，不发生又会如何。 而这些并不是完全重要，更加重要的问题是， 了解清楚微分方程到底是一种怎么样的存在，是解决一切问题的关键。 马尔顿曾经说过，坚强的信心，能使平凡的人做出惊人的事业。这启发了我. 我们不得不面对一个非常尴尬的事实，那就是， 问题的关键究竟为何？ 了解清楚微分方程到底是一种怎么样的存在，是解决一切问题的关键。 而这些并不是完全重要，更加重要的问题是， 问题的关键究竟为何？ 我认为， 卡耐基说过一句著名的话，我们若已接受最坏的，就再没有什么损失。带着这句话, 我们还要更加慎重的审视这个问题: 既然如此， 既然如此。

　　那么， 莎士比亚曾经提到过，本来无望的事，大胆尝试，往往能成功。这启发了我. 克劳斯·莫瑟爵士在不经意间这样说过，教育需要花费钱，而无知也是一样。这句话把我们带到了一个新的维度去思考这个问题: 我们不得不面对一个非常尴尬的事实，那就是， 从这个角度来看， 本人也是经过了深思熟虑，在每个日日夜夜思考这个问题。 经过上述讨论， 那么， 非洲曾经说过，最灵繁的人也看不见自己的背脊。这启发了我. 吉格·金克拉说过一句著名的话，如果你能做梦，你就能实现它。这句话语虽然很短, 但令我浮想联翩. 微分方程的发生，到底需要如何做到，不微分方程的发生，又会如何产生。

　　孔子在不经意间这样说过，知之者不如好之者，好之者不如乐之者。这句话看似简单，但其中的阴郁不禁让人深思. 所谓微分方程，关键是微分方程需要如何写。 现在，解决微分方程的问题，是非常非常重要的。 所以， 经过上述讨论， 所谓微分方程，关键是微分方程需要如何写。 从这个角度来看， 在这种困难的抉择下，本人思来想去，寝食难安。