第7.1讲 常微分方程与MATLAB求解

参数说明
（4）y0 包含初始条件的向量。 （5）options 用命令odeset设置的可选
积分参数. （6）p1,p2,… 传递给函数odefun的可选
参数。
[T,Y] = solver(odefun,tspan,y0)
在区间tspan=[t0,tf]上，从t0到tf，用 初始条件y0求解显式微分方程y’=f(t,y)。
y(xn1) yn1
则上式可记为
yn1 yn hf (xn , yn )
此即为求解初值问题的Euler方法，又称显式Euler方法。 27
龙格－库塔法
二阶龙格－库塔法
yn1 yn h(c1K1 c2 K 2 ) K1 f (xn , yn )
K 2 f (xn 2h, yn 21hK1 )
它每步求解yk+1需要解一个隐式方程。
24
欧拉（Euler）方法
欧拉（Euler）方法
在x= x0 处，用差商代替导数：
y(x0 )

y(x1) y(x0 ) x1 x0

y(x1) h
y(x0 )
由 y(x0 ) f (x0 , y0 ),
y(x0 ) y0
得 y(x1) y0 hf (x0 , y0 ) y1
相邻两个节点的间距 h xi1 xi 称为步长，步 长可以相等，也可以不等。假定h为定数，称为定步长，
这时节点可表示为
xi x0 ih, i 1,2,, n
20
离散化
数值解法需要把连续性的问题加以离散化，从 而求出离散节点的数值解。
对常微分方程数值解法的基本出发点就是离散化。 其数值解法的基本特点：采用“步进式”， 即求解过程顺着节点排列的次序一步一步地向前推进，
x1 o
本次课程的学习目标
理解常微分方程的涵义
知
识 应用MATLAB解决工程常微分问题 构建工程问题模型（微分方程部分）
能 讨论逻辑思维能力培养 力 讨论演示表达能力培养
常微分方程的理论基础
微分方程分类
包含自变量、未知函数及未知函数的导数或微分 的方程称为微分方程。
在微分方程中, 自变量的个数只有一个, 称为常微 分方程。
而y=[y;y(1);y(2);…,y(m-1)]， n与m可以不等
求解具体ODE的基本过程
（2）运用数学中的变量替换： yn=y(n-1),yn-1=y(n-2),…,y2=y‘,y1=y， 把高阶（大于2阶）的方程（组）写成一阶 微分方程组：
y


y1 y2
yn
常微分方程的符号解
(6)若没有给定输出参量，则在命令窗口显示 解列表。若该命令找不到解析解，则返回 一警告信息，同时返回一空的sym对象。 这时，用户可以用命令ode23或ode45求 解方程组的数值解。
例1
y a2 y

y(0) 1

y(
/
a)

0
>> dsolve('D2y = -a^2*y,y(0) = 1,Dy(pi/a) = 0','x')
数值计算方法与 MATLAB
包能胜 汕头大学工学院
第七讲 常微分方程与MATLAB
常微分方程的工程来由
常微分方程的工程来由
问题： 科学技术和工程中许多问题是用微分方程的形 式建立数学模型。因此微分方程的求解有很实 际的意义。
dV
1 D 2 ( x)
dx 4
V (0) 0
微分方程数值方法的基本思想
微分方程数值方法的基本思想
对常微分方程初值问题的数值解法，就是要算出精 确解y(x)在区间a,b上的一系列离散节点
a x0 x1 xn1 xn b
处的函数值 y(x0 ), y(x1),, y(xn )的近似值
y0 , y1 ,, yn
31
常微分方程的数值解
常微分方程数值解
Matlab专门用于求解常微分方程的函 数，主要采用Runge-Kutta方法：
ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb
龙格－库塔法
基于龙格－库塔法， MATLAB求常微分方程数值解的函数， 一般调用格式为：
自变量的个数为两个或两个以上的微分方程叫偏 微分方程。
微分方程中出现的未知函数最高阶导数的阶数称 为微分方程的阶数。
7
常微分方程
y f (x, y) a x b
(1)

y(a)

y0
y f (x, y, y) a x b
(2)y(a) y0 , y(a)
其中 c1, c2, 2, 21 待定。
29
二阶龙格－库塔法

1 c1 c2 0
1

2

c22

0
即
1 2

c2 21

0

c1 c2 1

c 2 2


1 2
c2 21

1 2
方程组解不唯一，可令c2=a 0 ，则 c1 = 1-a ， 2 = 21 =1/(2a)
21
描述这类算法，要求给出用已知信息 yi , yi1 , yi2 ,, y0 计算 yi1的递推公式。
建立这类递推公式的基本方法是在这些节点上用数值
积分、数值微分、泰勒展开等离散化方法，对初值问
题
y f (x, y)

y(
x0
)

y0
中的导数 y 进行不同的离散化处理。
22
将参数p1,p2,p3,..等传递给函数odefun，再 进行计算。若没有参数设置，则令 options=[]。
求解具体ODE的基本过程
求解具体ODE的基本过程
（1）根据问题所属学科中的规律、定律、 公式，用微分方程与初始条件进行描述。 F(y,y’,y’’,…,y(n),t) = 0 y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1



f1 (t, f2 (t,

fn (t,
y) y)
y)
y1 (0) y0
y0


y2
(
0)



y1



yn
(0)

yn

（3）根据（1）与（2）的结果，编写能计算
导数的M-函数文件odefile。
ans=
cos(a*x)
>> dsolve(‘D2y = -a^2*y’, ‘y(0) = 1,Dy(pi/a) = 0’)
例2 u v v u
>> [u,v] = dsolve('Du=v,Dv=u')
u= C1*exp(-t)+C2*exp(t) V= -C1*exp(-t)+C2*exp(t)
用参数options（用命令odeset生成）设置属 性（代替了缺省的积分参数），再进行操作。 常用的属性包括相对误差值RelTol（缺省值为 1e-3）与绝对误差向量AbsTol（缺省值为每 一元素为1e-6）。
[T,Y] = solver(odefun,tspan,y0,options,p1,p2…)
23
单步法和多步法
单步法：在计算yi+1 时只利用y i 多步法：在计算yi+1 时不仅利用y i , 还要利用 yi−1,
yi−2,…, k步法：在计算yi+1 时要用到yi,yi−1,…,yi−k+1 显式格式:可写成:yk+1=yk+hΦf(xk，yk;h） 隐式格式:yk+1=yk+hΦf（xk，yk,yk+1;h）
34
常微分方程数值解
[T,Y] = solver(odefun,tspan,y0) [T,Y] = solver(odefun,tspan,y0,options)
[T,Y] = solver(odefun,tspan,y0,options,p1,p2…)
参数说明
（1）solver为命令 ode45,ode23,ode113,ode15s, ode23s,ode23t,ode23tb之一。
对于标量t与列向量y，函数f=odefun(t,y) 必须返回一f(t,y)的列向量f。
解矩阵Y中的每一行对应于返回的时间列向量 T中的一个时间点。
要获得问题在其他指定时间点t0,t1,t2,…上 的解，则令tspan=[t0,t1,t2,…,tf]（要求是单 调的）。
[T,Y] = solver(odefun,tspan,y0,options)
[t,y]=ode23('fname',tspan,y0)
[t,y]=ode45('fname',tspan,y0)
其中fname是定义f(t,y)的函数文件名，该函数文件
必须返回一个列向量。tspan形式为[t0,tf],表示求
解区间。y0是初始状态列向量。t和y分别给出时间
向量和相应的状态向量。
（2）odefun 为常微分方程y’=f(x,y)， 或为包含一混合矩阵的方程 (x,y)*y’=f(x,y).
（3）tspan 积分区间（即求解区间）的向 量tspan=[t0,tf]。要获得问题在其他指定 时间点t0,t1,t2,…上的解，则令 tspan=[t0,t1,t2,…,tf] （要求是单调的）。
(1)，(2)式称为初值问题.