matlab解高阶微分方程程序

实验二： 微分方程模型Matlab求解与分析

一、实验目的

[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；

[2] 熟悉MATLAB软件关于微分方程求解的各种命令；

[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程；

[4] 熟悉离散 Logistic模型的求解与混沌的产生过程。

二、实验原理

1. 微分方程模型与MATLAB求解

解析解

用MATLAB命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程（组）的解析解。其中‘eqni'表示第i个微分方程，Dny表示y的n阶导数,默认的自变量为t。

（1） 微分方程

例1 求解一阶微分方程 21ydxdy

(1) 求通解

输入：

dsolve('Dy=1+y^2')

输出：

ans =

tan(t+C1)

（2）求特解

输入：

dsolve('Dy=1+y^2','y(0)=1','x')

指定初值为1，自变量为x

输出：

ans =

tan(x+1/4*pi)

例2 求解二阶微分方程 221()04(/2)2(/2)2/xyxyxyyy

原方程两边都除以2x，得211(1)04yyyxx

输入:

dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')

ans =

- (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) +

(exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2))

试试能不用用simplify函数化简

输入: simplify(ans)

ans =

2^(1/2)*pi^(1/2)/x^(1/2)*sin(x)

（2）微分方程组

例3 求解 df/dx=3f+4g; dg/dx=-4f+3g。

MATLAB中ODE的使⽤

ode23 解⾮刚性微分⽅程，低精度，使⽤Runge-Kutta法的⼆三阶算法。

ode45 解⾮刚性微分⽅程，中等精度，使⽤Runge-Kutta法的四五阶算法。

ode113 解⾮刚性微分⽅程，变精度变阶次Adams-Bashforth-Moulton PECE算法。

ode23t 解中等刚性微分⽅程，使⽤⾃由内插法的梯形法则。

ode15s 解刚性微分⽅程，使⽤可变阶次的数值微分（NDFs）算法。

ode23s 解刚性微分⽅程，低阶⽅法，使⽤修正的Rosenbrock公式。

ode23tb 解刚性微分⽅程，低阶⽅法，使⽤TR-BDF2⽅法，即Runger-Kutta公式的第⼀级采⽤梯形法则，第⼆

级采⽤Gear法。

[t，YY]=solver('F',tspan,Yo)

解算ODE初值问题的最简调⽤格式。

solver指上⾯的指令。

tspan=[0,30]; ％时域t的范围

y0=[1;0]; ％y（1）y(2)的初始值

[tt,yy]=ode45(@DyDt,tspan,y0);

plot(tt,yy(:,1)),title('x(t)')

function ydot=DyDt(t,y)

ydot=[y(2); 2*(1-y(1)^2)*y(2)-y(1)]

刚性⽅程：刚性是指其Jacobian矩阵的特征值相差⼗分悬殊。在解的性态上表现为，其中⼀些解变化缓慢，另

⼀些变化快，且相差较悬殊，这类⽅程常常称为刚性⽅程，⼜称为Stiff⽅程。

刚性⽅程和⾮刚性⽅程对解法中步长选择的要求不同。

刚性⽅程⼀般不适合由ode45这类函数求解，⽽应该采⽤ode15s等。

如果不能分辨是否是刚性⽅程，先试⽤ode45，再⽤ode15s。

[t，YY，Te，Ye,Ie] = solver('F',tspan,Yo,options,p1,p2,…)

解算ODE初值问题的最完整调⽤格式。

为了能够解出⽅程，要⽤指令odeset确定求解的条件和要求。在MATLAB中，求解⽅程组的指令都有默认的求

matlab求解微分方程数值解与解析解

微分方程是数学中的重要内容，它描述了物理、工程、经济等领域中的许多现象和问题。在实际应用中，我们经常需要求解微分方程的解析解或数值解。本文将以Matlab为工具，探讨如何求解微分方程并对比解析解与数值解的差异。

一、引言

微分方程是描述自然界中许多现象和问题的数学语言，它包含了未知函数及其导数与自变量之间的关系。微分方程的求解可以帮助我们了解问题的性质和变化规律，并为实际应用提供参考。在许多情况下，微分方程的解析解很难求得，这时我们可以利用计算机进行数值求解。

二、微分方程的数值解法

1.欧拉法

欧拉法是最简单的数值求解微分方程的方法之一。它通过将微分方程转化为差分方程，然后利用离散的点逼近连续的解。具体步骤如下：

（1）将微分方程转化为差分方程，即用近似的导数代替真实的导数；

（2）选择初始条件，即确定初始点的值；

（3）选择步长和求解区间，即确定求解的范围和步长； （4）使用迭代公式计算下一个点的值；

（5）重复步骤（4），直到达到指定的求解区间。

2.改进的欧拉法

欧拉法存在精度较低的问题，为了提高精度，可以使用改进的欧拉法。改进的欧拉法是通过使用两次导数的平均值来计算下一个点的值，从而提高了数值解的精度。

3.龙格-库塔法

龙格-库塔法是一种常用的数值求解微分方程的方法，它通过使用多个点的导数来逼近连续解。龙格-库塔法的步骤如下：

（1）选择初始条件和步长；

（2）使用迭代公式计算下一个点的值；

（3）计算下一个点的导数；

（4）根据导数的值和步长计算下一个点的值；

（5）重复步骤（3）和（4），直到达到指定的求解区间。

龙格-库塔法的精度较高，适用于求解一阶和高阶微分方程。

三、微分方程的解析解

解析解是指能够用公式或函数表示的方程的解。有些微分方程具有解析解，可以通过数学方法求得。例如，一阶线性常微分方程和某些特殊类型的二阶微分方程等。解析解的优势在于精确性和直观性，能够帮助我们深入理解问题的本质。

第五章 控制系统仿真

§ 5.2 微分方程求解方法

以一个自由振动系统实例为例进行讨论。如下图 1 所示弹簧 - 阻尼系统，参数如下：

M=5 kg, b=1 N.s/m, k=2 N/m, F=1N

x

b

M F

k

图 1 弹簧-阻尼系统

假设初始条件为： t0 0 时，将 m 拉向右方，忽

略小车的摩擦阻力， x(0) 0m x(0) 0m/ s 求系

统的响应。

) 用常微分方程的数值求解函数求解包括 ode45、ode23、ode113、ode15s、ode23s 等。

wffc1.m myfun1.m

一、常微分方程的数值求解函数 ode45 求解

解：系统方程为

m x b x kx F

这是一个单变量二阶常微分方程。

将上式写成一个一阶方程组的形式，这是函数

ode45 调用规定的格式。

令： x(1) x （位移）

x(2) x x(1) （速度）

上式可表示成：

x(1) x(2) x(2)

x(2) x 1 10* x(2) 20* x(1)

下面就可以进行程序的编制。

%写出函数文件 myfun1.m

function xdot=myfun1(t,x)

xdot=[x(2);1-10*x(2)-20*x(1)];

% 主程序 wffc1.m t=[0 30];

x0=[0;0];

[tt,xx]=ode45(@myfun1,t,x0);

plot(tt,yy(:,1),':b',tt,yy(:,2),'-r')

legend(' 位移 ',' 速度 ')

title(' 微分方程的解 x(t)')

二、方法 2：

m x b x kx F

X(s) 1

G(s)

5s2 s 2 F (s)

%用传递函数编程求解 ksys1.m

num=1;

den=[5 1 2];

%printsys(num,den)

%t=0:0.1:10;