# Matlab求解常微分方程

# 理论知识 ## 概念 微分方程：含导数或微分的方程。
解：满足微分方程的函数。
特解/通解：特解指的是满足微分方程的某一个解；通解指的是满足微分方程的一组解。 阶：微分方程中导数或微分的最高阶数。
线性/非线性：（几何意义：叠加原理）方程中的函数和它的各阶导数都是一次方为线性微分方程，否则为非线性。例： - y'=sin(x)*y 线性 - y'=y^2 非线性
齐次/非齐次：（代数意义：次数）齐次微分方程中不含常数项,也不含仅由x的各种运算组合构成的项（比如4xx,sinx等），否则为非齐次。
常微分/偏微分：未知函数是一元函数的，叫常微分方程；未知函数是多元函数的叫做偏微分方程。
初值问题和边值问题：附加条件中未知函数及其导数的独立变量取值相同，则为初值问题；附加条件中未知函数及其导数的独立变量取值不同，则为边界值问题。例： - y(0)=1,y'(0)=2 初值问题 - y(0)=1,y'(1)=2 边值问题 ## 解法 在Matlab中解微分方程，方法有两种： - 解析解方法：严格按照公式逻辑推导得到的，具有基本的函数形式。 - 数值解方法：采用某种计算方法，在特定的条件下得到的一个近似数值结果，如有限元法，数值逼近法，插值法等等。
解析解方法可回顾《高等数学》中相关公式。 数值解方法可回顾《数值分析》中相关解法。

Matlab微分方程求解

解析解

解析解的存在

由Abel-Ruffini定理，四次及以下的多项式代数方程是能求出根的解析解的，即低阶常系数线性微分方程有一般意义下的解析解。

非线性微分方程只能用数值解法求解，即使看起来很简单的非线性微分方程也是没有解析解的，只有极特殊的非线性微分方程解析可解。

解析解的解法

利用dsolve函数

S = dsolve(eqn) S = dsolve(eqn,cond) S = dsolve(,Name,Value) [y1,...,yN] = dsolve()

实例

例1 输入信号为u(t)=exp(-5t)cos(2t+1)+5，求微分方程diff(y,4)+10diff(y,3)+35diff(y,2)+50diff(y)+24y=5diff(u,t,2)+4diff(u,t)+2u的通解。初值条件：y(0)=3,y1(0)=2,y2(0)=0,y3(0)=0，求方程的特解。（约定y1代表一阶导数，以此类推）

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
% y=dsolve(f1,f2,...,fm) 默认自变量为t
% y=dsolve(f1,f2,...,fm,'x') 指明自变量
% f可由字符串表示也可由符号表达式表示
 
%用字符串表达式，新版本会被移除
syms t; 
u=exp(-5*t)*cos(2*t+1)+5;
uu=5*diff(u,t,2)+4*diff(u,t)+2*u;
y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)]); %求解
y=simplify(y) %化简
 
%用符号表达式，推荐
syms y(t); 
eqn=diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==uu;
y=dsolve(eqn);
y=simplify(y) 
 
%检验
diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y-uu
 
%求特解并绘图
syms y(t); 
eqn=diff(y,4)+10*diff(y,3)+35*diff(y,2)+50*diff(y)+24*y==uu;
y1=diff(y);y2=diff(y,2);y3=diff(y,3);y4=diff(y,4); 
% 需要引入中间变量
% 用字符串求解的情况，不需要引入中间变量。'y(0)=3','Dy(0)=2','D2y(0)=0'...
cond=[y(0)==3,y1(0)==2,y2(0)==0,y3(0)==0];
z=dsolve(eqn,cond);
z=simplify(z)
ezplot(z,[0,5])
fplot(z,[0,5])
% ezplot(fun,[xmin,xmax])绘制fun(x)在以下域上的图形:xmin<x<xmax
% ezplot适用隐式函数
% 推荐fplot代替ezplot
double(subs(z,t,5)) 
% 验证图像，对自变量赋值，并求小数解
 
% 改变初值条件，再求特解，用字符串表示方式
z1=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',char(uu)],...
    'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5') %求解
fplot(z1,[0,5])

运行结果：

y = C1exp(-4t) - (547exp(-5t)sin(2t + 1))/520 - (343exp(-5t)cos(2t + 1))/520 + C2exp(-3t) + C3exp(-2t) + C4*exp(-t) + 5/12

ans =0

z = 19exp(-t) - (69exp(-2t))/2 + (73exp(-3t))/3 - (25exp(-4t))/4 + (97exp(-t)sin(1))/60 - (51exp(-2t)sin(1))/13 + (5exp(-3t)*sin(1))/8......

ans = 0.5679

z1 = (exp(-5t)(445cos(2t + 1) - 65exp(5t) + 102sin(2t + 1)))/26 - (exp(-5t)(537cos(2t + 1) - 40exp(5t) + 15sin(2t + 1)))/24 - (exp(-5t)(25exp(5t) - 542cos(2t + 1) ......

例2 求解有复数极点的微分方程diff(y,5)+5diff(y,4)+12diff(y,3)+16diff(y,2)+12diff(y)+4y=diff(u)+3u，输入信号：u(t)=sin(t)，初值条件：y(0)=0,y1(0)=0,y2(0)=0,y3(0)=0,y4(0)=0

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
% 有复数极点的微分方程
syms y(t) u(t); 
u=sin(t);
uu=diff(u)+3*u;
eqn=diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y==uu;
y1=diff(y);y2=diff(y,2);y3=diff(y,3);y4=diff(y,4); % 需要引入中间变量
cond=[y(0)==0,y1(0)==0,y2(0)==0,y3(0)==0,y4(0)==0];
y=dsolve(eqn,cond);
y=simplify(y)
fplot(y,[0,30])
% 检验
diff(y,5)+5*diff(y,4)+12*diff(y,3)+16*diff(y,2)+12*diff(y)+4*y-uu

运行结果：

y = exp(-t) - cos(t)/5 - (2sin(t))/5 - (4exp(-t)cos(t))/5 + (11exp(-t)sin(t))/10 - (texp(-t)*cos(t))/2

ans =0

数值解

概述

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

x ' ( t ) = f ( t , x ( t ) )

x(t)=[x1(t);x2(t);...xn(t)] 状态向量 f(.)=[f1(.);f2(.);...fn(.)] 任意非线性函数 x0=[x1(t0);x2(t0);...xn(t0)] 初始状态 [t0,tn] 时间区间（tn为终止时间） ### 优化措施 选择适当的步长h：不能太大和太小，太大精度不够，太小减慢计算精度，增加计算次数从而增加累积误差。 改进近似算法精度：Euler算法将原始积分进行梯形近似，精度低；精度较高的是Runge-Kutta法和Adams法。 采用变步长方法：采用自适应步长算法。 ### 解法 1、手动编写算法。例如：

四阶定步长Runge-Kutta算法matlab代码如下：

2、ode系列函数：最常用ode45() 介绍

ode45()实现了变步长四阶五级 Runge-Kutta-Felhberg算法，可以使用变步长的方法求解微分方程。ode系列函数功能总结如下：

用法

函数调用：

[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(___) 微分方程/组描述：

function xd=funname(t,x) % 不需要附加参数格式 function xd=funname(t,x,p1,p2,...,pm) % 可以使用附加参数 控制条件：使用odesrt()函数创建或修改 options 结构体。

options=odeset('RelTol',1e-7) options=odeset;options.RelTol=1e-7 例：options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) 指定 1e-5 的相对误差容限、打开求解器统计信息的显示，并指定输出函数 @odeplot 在计算时绘制解。

检验

修改仿真控制参数，例如：将‘RelTol’或'AbsTol'选项设置成一个更小的值，观察结果是否和上次一致，如有不可接受的差异，应考虑进一步减小误差值。 选择不同的微分方程求解算法。