次课数值微积分与常微分方程求解-修订.pptVIP

数值微积分与常微分方程求解 【学习目标】 ● 掌握微分与积分的数值计算方法。 ● 掌握常微分方程的数值求解方法。 数 值 微 分 数值差分与差商 数值微分的实现 1．基本思想是先用逼近或拟合等方法将已知数据在一定范围内的近似函数求出，再对此近似函数进行微分。 比如用多项式函数g(x)对f(x)进行逼近，然后用g(x)在点x处的导数作为f(x)在点x处的导数。该种方法一般只用在低阶数值微分。 2．用f(x)在点x处的某种差商作为其导数。 在MATLAB中，提供计算向前差分的函数diff。 ● DX?=?diff(X)：计算向量X的向前差分，DX(i)?=?X(i+1) ?X(i)，i?=?1,2, …,n??1。 ● DX?=?diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)?=?diff(diff(X))。 ● DX?=?diff(A,n,dim)：计算矩阵A的n阶差分，dim?=?1时（默认状态），按列计算差分；dim?=?2，按行计算差分。 函数diff计算的是向量元素间的差分，故所得输出比原向量少了一个元素。 【例】生成一个5阶魔方矩阵，按列进行差分运算。 M=magic(5) M= 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 DM=diff(M) DM= 6 ?19 6 6 1 ?19 1 6 6 6 6 6 6 1 ?19 1 6 6 ?19 6 可以看出，diff函数对矩阵的每一列都进行差分运算，因而结果矩阵的列数是不变的，只有行数减1。矩阵DM第3列值相同，表明原矩阵第3列是等间距的。 【例】设f(x)?=?sinx，用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f?(x)的图像。 为确定计算数值导数的点，假设在[0， π]区间内以π/24为步长求数值导数。下面用3种方法求f(x)在这些点的导数。x=0:pi/24:pi; %用5次多项式p拟合f(x),并对拟合多项式p求导数dp在假设点的函数值 p=polyfit(x,sin(x),5); dp=polyder(p); dpx=polyval(dp,x); dx=diff(sin([x,pi+pi/24]))/(pi/24); gx=cos(x); plot(x,dpx,x,dx,o,x,gx,+); 定积分的数值求解实现 1．自适应辛普生法 quad函数和quadl函数 函数的调用格式为 [I,n]=quad(@fname,a,b,tol,trace) [I,n]=quadl(@fname,a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，默认时取tol?=?10-6。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace?=?0。返回参数I即定积分值，n为被积函数的调用次数。 一般情况下，quadl函数调用的步数明显小于quad函数，而且精度更高，从而保证能以更高的效率求出所需的定积分值。 【例】分别用quad函数和quadl函数求椭圆积分 的近似值，并在相同的积分精度下，比较函数的调用次数。 调用函数quad求定积分： format long; fx=inline(1./sqrt(1+X.^4)); %定义一个语句函数 [I,n]=quad(fx,0,1,1e-10) %注意函数名不加@号 I= 0.927037338654481 n= 113 调用函数quadl求定积分： [I,n]=quadl(fx,0,1,1e-10) I= 0.927037338650659 n= 48 2．高斯-克朗罗德法 该函数的调用格式为 [I,err]?=?quadgk(@fname,a,b) 其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是?Inf或Inf，也可以是复数。如果积分上下限是复数，则quadgk在复平面上求积分。 【例】求定积分 （1）建立被积函数文件feln.m。 functio