matlab常微分方程数值解法精品PPT课件
合集下载
微分方程的Matlab求解ppt课件
取t0=0，tf=10，建立程序如下： t0=0;tf=10; [t,y]=ode45('eq3',[t0 tf],[0 0]); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,'-') hold on plot(y(:,1),y(:,2),'*')
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
例4
d2x dt 2
1000(1
x2
)
dx dt
x
0
x(0) 2; x'(0) 0
解: 令 y1=x，y2=y1’
由待解 方程写 成的m-
ts=[t0，tf]， t0、tf为自
函数的 初值
变量的初
文件名 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
z = (-C1+exp(4*t)*C1-2*exp(4*t)+C2+exp(4*t)*C3)*exp(-2*t)
用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
例4
d2x dt 2
1000(1
x2
)
dx dt
x
0
x(0) 2; x'(0) 0
解: 令 y1=x，y2=y1’
由待解 方程写 成的m-
ts=[t0，tf]， t0、tf为自
函数的 初值
变量的初
文件名 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
z = (-C1+exp(4*t)*C1-2*exp(4*t)+C2+exp(4*t)*C3)*exp(-2*t)
用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s
matlab常微分方程 ppt课件
odeplot 随计算过程，显示解向
量
属性名 OutputSel
Refine
2020/4/10
参数设置
取值
含义
若不使用缺省设置，则
有效值： 正整数向 量 缺省值：[]
OutputFcn所表现的是那 些正整数指定的解向量 中的分量的曲线或数据 。若为缺省值时，则缺 省地按上面情形进行操
作
有效值： 正整数k>1 缺省值：k
y1 (0 ) y 0
y0
y
2
(
0
)
y1
yn
(
0
)
yn
• （3）根据（1）与（2）的结果，编写能计 算
• 导数的M-函数文件odefile。
• （4）将文件odefile与初始条件传递给求解
• 器Solver中的一个，运行后就可得到ODE
• 的、在指定时间区间上的解列向量y（其中包 2含020/4y/10及不同阶的导数）。
参数设置
属性名 NormControl Events
取值
含义
有效值： on、off 缺省值：
off
为‘on’时，控制解向量 范数的相对误差，使每 步计算中，满足： norm(e)<=max(RelTol*n orm(y),AbsTol)
有效值： 为‘on’时，返回相应的 on、off 事件记录
2020/4/10
例3
x' x 2
x(0)
1
创建函数function2, 保存在function2.m中
function f=function2(t,x)
f=-x.^2;
在命令窗口中执行 >> [t,x]=ode45('function2',[0,1],1); >> plot(t,x,'-',t,x,'o'); >> xlabel('time t0=0,tt=1'); >> ylabel('x values x(0)=1');
如何使用MATLAB求解微分方程(组)ppt课件
差,输出参数,事件等，可缺省。 9
使用ODE?时如何编 写微分方程 ？
方式一：带额外参数，使用时需对参数进行赋值
function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变 量
xdot=zeros(2,1);
%表明其为列向量
xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);
2
Where ?
工程控制
ODE
医学生理
航空航天
金融分析
3
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
4
When ?
当对问题进行建模后，有常微分方程 需要求解时。
在生物建模中，经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
5
方法 ode23
ode45
数 ode113
当无法藉由微积分技巧求 得解析解时，这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下，常 微分方程往往只能求解出其
数值解。
在数学中，刚性方程是指一 个微分方程，其数值分析的解 只有在时间间隔很小时才会稳 定，只要时间间隔略大，其解 就会不稳定。
目前很难去精确地去定义哪 些微分方程是刚性方程，但是 大体的想法是：这个方程的解
y(1)=x(2);
y1
y2
y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
1000
求解程序ห้องสมุดไป่ตู้键步骤
[t,y]=ode45('odefun',[3.9 4.0],[2 8])
y
《微分方程的数值解法maab四阶龙格—库塔法》PPT模板课件
k1 f (tn , yn )
k2
f (tn
1 2
h,
yn
h 2 k1)
k3
f (tn
1 2
h,
yn
h 2
k2
)
k 4 f (tn h , y n hk 3 )
四 阶 Runge-Kutta 法计算流程图
开始
h 初始条件：t
迭代次数：
；y
0
N
0
积分步长：
tn t0
for i = 1 : N
解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P 8 k 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
各种solver 解算指令的特点
解法指令 解题类 型
特点
ode45 非刚性 采用4、5阶Runge－Kutta法
适合场合 大多数场合的首选算法
ode23 非刚性 采用Adams算法
较低精度（10－3）场合
ode113 非刚性
ode23t ode15s
适度刚 性
刚性
ode23s 刚性 ode23tb 刚性
y2
(0)
(2.2) (2.3)
例：著名的Van der Pol方程
y (y21)y y0
令
y1y,y2y
Y
y y
1 2
降为一阶 Yyy12(y12y12)y2y1
初始条件
Y0
y1(0) y2(0)
y10 y20
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
科学计算与MATLAB
2012
第七讲 常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
作如下近似
yn y(xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn1 y(xn h)
y(xn )
hy(xn )
1 2
h2
y(xn )
y(xn )
x2 …. xn ….
y(x0) y(x1) y(x2) …. y(xn) ….
y0
y1 y2 …. yn ….
在xn节点上，微分方程可以写为
y(xn1) y(xn ) f y(xn ) , xn h
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , xn ) h
hf
y(xn ),
xn
1 2
h2 y(xn )
作如下近似
yn y(xn )
yn1 yn f yn , xn h
局部截断误差
y
y0
dy dx
x0 x x0
y0 f ( y0 , x0 ) x x0
此切线与x=x1交点纵坐标为：
y1 y0 f ( y0 , x0 ) x1 x0
=y0 f ( y0 , x0 ) h
y1
Qy((tt))
t0 t1 t2 t3 t4 t5 t6 t
从（x1,y1）作曲线y(x)在 (x1,y(x1))的切线的平行线，
作曲线y(x)在(xn,y(xn))的切
Q1
线的平行线，得直线方程：
Qy((tt))
y yn f ( yn , xn ) x xn
与x=xn+1交点纵坐标为：
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
yn1 yn f ( yn , xn ) h
折线近似曲线，n增大，误差变大
Qy
得x1
y1 f ( y1, x1) x x1
与x=x2交点纵坐标为：
y2 y1 f ( y1, x1) h
y1
Qy((tt))
y2
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
y1 y(x1)
按照相似的方法，从（xn,yn）Qy
yn y(xn )
2、欧拉近似方法
2.1 简单欧拉(L.Euler, 1707-1783)方法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。
y0 y1 y2 yn
欧拉数值算法递推公式构造 2.1.1 差分法
f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值
end
outy=y;
outx=x; %plot(x,y)%画出方程解的函数图
例题： y ' sin x y y(x0 ) 1, x0 0
2.1.2 折线法（几何意义） Qy 从（x0,y0）作曲线y(x)的切 线，得切线方程：
PointNum=100;
end if nargin<4 %y0默认值为0
y0=0;
end h=(xt-x0)/PointNum;%计算步长h x=x0+[0:PointNum]'*h;%自变量数组 y(1,:) = y0(:)';%将输入存为行向量，输入为列向量形式
for k = 1:PointNum f=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值
同样，在[x0,x2] ，积分采用矩形近似，得：
y(x2 ) y0
x2 f y(x), x dx
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
2.1.3 数值积分
对微分方程作从 x0 到 x 积分得：
y(x) dy y0
x
x0 f y
x
, xdx
y x y0
x x0
f y x, xdx
在[x0,x1] ，积分采用矩形近似，得：
y(x1) y0
x1 f y(x), x dx
x0
y0 f y(x0 ), x0 h
dt 2
简单问题可以求得解析解，多数实际问题靠数值求解。
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 yn，用以近似y(xn)
差分法就是用差商近似代替微商，即
y dy x dx
代入微分方程得到：
y(x x) y(x) f ( y(x), x) x
y(x x) y(x) f ( y(x), x)x
对于等间隔节点
x xn1 xn h xn1 xn h
n=0,1,2
可以得到：
xn y精确值 y近似值
x0 x1
具体求解过程为：
y1 y0 f ( y0 , x0 ) h y2 y1 f ( y1 , x1) h y3 y2 f ( y2 , x2 ) h
简单欧拉方法程序
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) %MyEuler 用前向差分的欧拉方法解微分方程 %fun 表示f（x，y） %x0,xt表示自变量的初值和终值 %y0表示函数在x0处的值,其可以为向量形式 %PointNum表示自变量在[x0,xt]上取的点数 if nargin<5 | PointNum<=0 %如果函数仅输入4个参数值，则PointNum默认值为100
2012
第七讲 常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
作如下近似
yn y(xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn1 y(xn h)
y(xn )
hy(xn )
1 2
h2
y(xn )
y(xn )
x2 …. xn ….
y(x0) y(x1) y(x2) …. y(xn) ….
y0
y1 y2 …. yn ….
在xn节点上，微分方程可以写为
y(xn1) y(xn ) f y(xn ) , xn h
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , xn ) h
hf
y(xn ),
xn
1 2
h2 y(xn )
作如下近似
yn y(xn )
yn1 yn f yn , xn h
局部截断误差
y
y0
dy dx
x0 x x0
y0 f ( y0 , x0 ) x x0
此切线与x=x1交点纵坐标为：
y1 y0 f ( y0 , x0 ) x1 x0
=y0 f ( y0 , x0 ) h
y1
Qy((tt))
t0 t1 t2 t3 t4 t5 t6 t
从（x1,y1）作曲线y(x)在 (x1,y(x1))的切线的平行线，
作曲线y(x)在(xn,y(xn))的切
Q1
线的平行线，得直线方程：
Qy((tt))
y yn f ( yn , xn ) x xn
与x=xn+1交点纵坐标为：
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
yn1 yn f ( yn , xn ) h
折线近似曲线，n增大，误差变大
Qy
得x1
y1 f ( y1, x1) x x1
与x=x2交点纵坐标为：
y2 y1 f ( y1, x1) h
y1
Qy((tt))
y2
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
y1 y(x1)
按照相似的方法，从（xn,yn）Qy
yn y(xn )
2、欧拉近似方法
2.1 简单欧拉(L.Euler, 1707-1783)方法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。
y0 y1 y2 yn
欧拉数值算法递推公式构造 2.1.1 差分法
f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值
end
outy=y;
outx=x; %plot(x,y)%画出方程解的函数图
例题： y ' sin x y y(x0 ) 1, x0 0
2.1.2 折线法（几何意义） Qy 从（x0,y0）作曲线y(x)的切 线，得切线方程：
PointNum=100;
end if nargin<4 %y0默认值为0
y0=0;
end h=(xt-x0)/PointNum;%计算步长h x=x0+[0:PointNum]'*h;%自变量数组 y(1,:) = y0(:)';%将输入存为行向量，输入为列向量形式
for k = 1:PointNum f=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值
同样，在[x0,x2] ，积分采用矩形近似，得：
y(x2 ) y0
x2 f y(x), x dx
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
2.1.3 数值积分
对微分方程作从 x0 到 x 积分得：
y(x) dy y0
x
x0 f y
x
, xdx
y x y0
x x0
f y x, xdx
在[x0,x1] ，积分采用矩形近似，得：
y(x1) y0
x1 f y(x), x dx
x0
y0 f y(x0 ), x0 h
dt 2
简单问题可以求得解析解，多数实际问题靠数值求解。
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 yn，用以近似y(xn)
差分法就是用差商近似代替微商，即
y dy x dx
代入微分方程得到：
y(x x) y(x) f ( y(x), x) x
y(x x) y(x) f ( y(x), x)x
对于等间隔节点
x xn1 xn h xn1 xn h
n=0,1,2
可以得到：
xn y精确值 y近似值
x0 x1
具体求解过程为：
y1 y0 f ( y0 , x0 ) h y2 y1 f ( y1 , x1) h y3 y2 f ( y2 , x2 ) h
简单欧拉方法程序
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) %MyEuler 用前向差分的欧拉方法解微分方程 %fun 表示f（x，y） %x0,xt表示自变量的初值和终值 %y0表示函数在x0处的值,其可以为向量形式 %PointNum表示自变量在[x0,xt]上取的点数 if nargin<5 | PointNum<=0 %如果函数仅输入4个参数值，则PointNum默认值为100