第5章 数值微积分与常微分方程求解
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值微积分与常微分方程求解
数值微分
diff函数：用f(x)在点x处得某种差商做为其导数。 在MATLAB中，没有直接提供求数值导数的 函数，
只有计算向前差分的函数diff，调用格式
• •
DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i) DX=diff(X)：计算向量X的n阶向前差分，例如， diff(X,2)=diff(diff(X)) DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时（默认
例：用trapz函数计算定积分
2 .5
1
1 dx 2 1 x
多重积分的数值求解：二重积分和三重积分函
数的调用格式为
dblquad (fun,a,b,c,d,tol)
triplequad (fun,a,b,c,d,e,f,tol)
其中，fun为被积函数，[a,b]为x的积分区域， [c,d]为y的积分区域，[e,f]为z的积分区域，参 数tol，trace的用法与函数quad完全相同
•
状态），按列计算差分；dim=2时，按行计算差分
注：对于向量的微分，函数diff计算的是向量元素间的差分
例：设f(x)=sin(x)，用不同的方法求函数f(x)的数值
导数，并在同一坐标系中作出f’(x)的图像
解：为确定计算数值导数的点，假设在[0,20]区间
内以pi为步长求数值导数。两种方法分别为
求解函数 ode23 oder45 oder113 ode23t
采用方法 2-3阶Runge-Kutta算法，低精度 4-5阶Runge-Kutta算法，中精度 Adam算法，精度可到10-3~10-6 梯形算法
使用场合 非刚性 非刚性 非刚性 适度刚性
Ode15s
Ode23s ode23tb
Gear’s反向数值微分算法，中精度 刚性
[I, err]=quadgk[@fname,a,b]
其中，err返回近似误差范围，其他参数的含义与
用法与quad函数相同，积分上下限可以是-Inf或
Inf，也可以是复数（在复平面上求积分）
梯形积分法：对由表格形式定义的Leabharlann Baidu数关系的求
定积分问题用梯形积分函数trapz，调用格式为 trapz(Y)：若Y是向量，则从1开始取单位步长， 以Y为函数值计算积分值；若Y是一矩阵，则计算 Y的每一列的积分
• •
用diff函数直接求f(x)在假设点的数值导数 先求出导函数f’(x)=cos(x)，然后直接求f’(x)在假设点的导数
数值积分
利用牛顿-莱布尼兹公式可以精确地计算定积分的 值，但仅适用于被积函数的原函数能用初等函数 表达出来的情形，大多数实际问题的积分需要用 数值方法求出近似结果
定积分的数值求解：在MATLAB中可以用quad或quad1来 进行数值积分 自适应simpson法，调用格式 [I, n]=quad(@fname, a,b,tol,trace) [I, n]=quadl(@fname, a,b,tol,trace) 其中，fname是被积函数名；a,b是定积分的上限和下限； tol用来控制积分精度，默认是取10-6; trace控制是否展现 积分过程，若取非0则展现，若取0则不展现，默认时取0； 返回参数I即定积分值；n为被积函数的调用次数
其中t和y分别给出时间向量和相应的状态向量；solver为求
常微分方程数值解的函数（下表函数之一）；fname是定义
f(t,y)的函数文件名，该函数文件必须返回一个列向量； tspan形式为[t0,tf]表示求解区间；y0是初始状态列向量；
options设置求解属性，常用的属性包括相对误差
值’RelTol’(默认为10-3)与绝对误差向量’AbsTol’(默认为 10-6)
2阶Rosebrock算法，低精度 梯形算法，低精度 刚性 刚性
例：设有如下初值问题，试求其数值解，并与精确
p 解相比较，精确解为( x) = ( x + ) / cos x y 2 p ¢= ytgx + sec x, 0 ＃x 1, y y = x =0 2
(1)
建立函数文件
(2)
求解微分方程
例：某非刚性物体的运动方程为
x x yz y ( y z ) z xy y z
初始条件为x(0)=0, y(0)=0, z(0)=ε;取β =8/3, ρ =28, σ =10, 试绘制系统相平面图
(1)
建立模型的函数文件
例：计算二重积分
1
2
1 2
e
x2 / 2
sin(x y)dxdy
2
例：计算三重积分
0 0
1
0
4 xze
z 2 y x2
dxdydz
常微分方程的数值求解
MATLAB提供了多个求解常微分方程数值解的函数， 一般调用格式为
[t,y]=solver(fname,tspan,y0[,options])
例：求
S
0
x sin(x) dx 1 cos(x)
(1)建立被积函数
(2)调用数值积分函数quad求定积分
例：分别用quad函数和quadl函数求椭圆如下积 分的近似值，并在相同的积分精度下，比较两个 函数的调用次数
1
1 1 x
4
0
dx
Gauss-Kronrod法：用来求震荡函数的定积分，调 用格式
(2)
解微分方程 绘制系统相平面图
(3)
数值微分
diff函数：用f(x)在点x处得某种差商做为其导数。 在MATLAB中，没有直接提供求数值导数的 函数，
只有计算向前差分的函数diff，调用格式
• •
DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i) DX=diff(X)：计算向量X的n阶向前差分，例如， diff(X,2)=diff(diff(X)) DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时（默认
例：用trapz函数计算定积分
2 .5
1
1 dx 2 1 x
多重积分的数值求解：二重积分和三重积分函
数的调用格式为
dblquad (fun,a,b,c,d,tol)
triplequad (fun,a,b,c,d,e,f,tol)
其中，fun为被积函数，[a,b]为x的积分区域， [c,d]为y的积分区域，[e,f]为z的积分区域，参 数tol，trace的用法与函数quad完全相同
•
状态），按列计算差分；dim=2时，按行计算差分
注：对于向量的微分，函数diff计算的是向量元素间的差分
例：设f(x)=sin(x)，用不同的方法求函数f(x)的数值
导数，并在同一坐标系中作出f’(x)的图像
解：为确定计算数值导数的点，假设在[0,20]区间
内以pi为步长求数值导数。两种方法分别为
求解函数 ode23 oder45 oder113 ode23t
采用方法 2-3阶Runge-Kutta算法，低精度 4-5阶Runge-Kutta算法，中精度 Adam算法，精度可到10-3~10-6 梯形算法
使用场合 非刚性 非刚性 非刚性 适度刚性
Ode15s
Ode23s ode23tb
Gear’s反向数值微分算法，中精度 刚性
[I, err]=quadgk[@fname,a,b]
其中，err返回近似误差范围，其他参数的含义与
用法与quad函数相同，积分上下限可以是-Inf或
Inf，也可以是复数（在复平面上求积分）
梯形积分法：对由表格形式定义的Leabharlann Baidu数关系的求
定积分问题用梯形积分函数trapz，调用格式为 trapz(Y)：若Y是向量，则从1开始取单位步长， 以Y为函数值计算积分值；若Y是一矩阵，则计算 Y的每一列的积分
• •
用diff函数直接求f(x)在假设点的数值导数 先求出导函数f’(x)=cos(x)，然后直接求f’(x)在假设点的导数
数值积分
利用牛顿-莱布尼兹公式可以精确地计算定积分的 值，但仅适用于被积函数的原函数能用初等函数 表达出来的情形，大多数实际问题的积分需要用 数值方法求出近似结果
定积分的数值求解：在MATLAB中可以用quad或quad1来 进行数值积分 自适应simpson法，调用格式 [I, n]=quad(@fname, a,b,tol,trace) [I, n]=quadl(@fname, a,b,tol,trace) 其中，fname是被积函数名；a,b是定积分的上限和下限； tol用来控制积分精度，默认是取10-6; trace控制是否展现 积分过程，若取非0则展现，若取0则不展现，默认时取0； 返回参数I即定积分值；n为被积函数的调用次数
其中t和y分别给出时间向量和相应的状态向量；solver为求
常微分方程数值解的函数（下表函数之一）；fname是定义
f(t,y)的函数文件名，该函数文件必须返回一个列向量； tspan形式为[t0,tf]表示求解区间；y0是初始状态列向量；
options设置求解属性，常用的属性包括相对误差
值’RelTol’(默认为10-3)与绝对误差向量’AbsTol’(默认为 10-6)
2阶Rosebrock算法，低精度 梯形算法，低精度 刚性 刚性
例：设有如下初值问题，试求其数值解，并与精确
p 解相比较，精确解为( x) = ( x + ) / cos x y 2 p ¢= ytgx + sec x, 0 ＃x 1, y y = x =0 2
(1)
建立函数文件
(2)
求解微分方程
例：某非刚性物体的运动方程为
x x yz y ( y z ) z xy y z
初始条件为x(0)=0, y(0)=0, z(0)=ε;取β =8/3, ρ =28, σ =10, 试绘制系统相平面图
(1)
建立模型的函数文件
例：计算二重积分
1
2
1 2
e
x2 / 2
sin(x y)dxdy
2
例：计算三重积分
0 0
1
0
4 xze
z 2 y x2
dxdydz
常微分方程的数值求解
MATLAB提供了多个求解常微分方程数值解的函数， 一般调用格式为
[t,y]=solver(fname,tspan,y0[,options])
例：求
S
0
x sin(x) dx 1 cos(x)
(1)建立被积函数
(2)调用数值积分函数quad求定积分
例：分别用quad函数和quadl函数求椭圆如下积 分的近似值，并在相同的积分精度下，比较两个 函数的调用次数
1
1 1 x
4
0
dx
Gauss-Kronrod法：用来求震荡函数的定积分，调 用格式
(2)
解微分方程 绘制系统相平面图
(3)