第七章微积分数值计算.ppt

X=[1 3 5 2 4]; Y= gradient(X) Y = 2.0000 2.0000 -0.5000 -0.5000 2.0000 Y=gradient(X,2) Y= 1.0000 1.0000 -0.2500 -0.2500 1.0000 即 两边用前向和后向差分，中间用中心差分 7.1.4 拉普拉斯算子4*del2 由于内部算法的原因： U=4*del2(v,h)，对1维向量v以步长h求拉普拉斯算子时，返回一相同维数的向量U，且 默认的步长为1。 U=4*del2(v,h1,h2)，对矩阵v，横向(x方向)以步长h1，纵向(y方向)以步长 h2计算拉普拉斯算子。 4*del2(U) ans = 4 4 4 4 4 4 4 4 4 4 4 4 [x,y]= meshgrid(1:4,1:3); U = x.*x+y.*y U = 2 5 10 17 5 8 13 20 10 13 18 25 首先编写vdp.m如下： function fy=vdp(x,y) %注意表示一个微分方程组的函数必须有自变量和因变量两个输入变量(向量) 。 fy=[y(2);7*(1-y(1).^2).*y(2)-y(1)]; 计算微方程组 y0=[1;0]； [x,y]=ode45(@vdp,[0,40],y0); % x 为程序在计算过程中所取的步长，y(:,1)=y，y(:,2)=y’，y(:,3)=y”…… y=y(:,1);dy=x(:,2); plot(x,y,x,dy,’r’) 7.1 数值微分 7.2 数值积分 7.3 常微分方程的数值解法 7.1 数值微分 实际问题常需计算函数的导数或积分值。但很多情况下，函数关系难以准确表示；即使能使用解析式准确表示，表示式却很复杂，不能用于实际计算。本章介绍数值计算导数或积分的实用方法。 7.1.1 差分和差商 根据导数的定义 其中，?x和?y分别称为自变量x和因变量y的增量，也称之为差分。可以用差分的商 作微商（导数）的近似。数值微分就是用函数值的线性组合近似函数在某点的导数值。自变量x的步长一般取定值。 首先在xi处对函数进行泰勒展开， 根据不同的组合方式可以得到精度不同的差分公式。以函数的一阶导数为例 ： 微分公式 表达式 截断误差 两点前向 O(?x) 两点后向 O(?x) 三点中心 O(?x2) 三点前向 O(?x2) 三点后向 O(?x2) 五点中心 O(?x4) …… …………………… …… 精度为O(?X2)的高阶中心差分算法 精度为O(?X4)的高阶中心差分算法 7.1.2 数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff和梯度函数gradient。 diff调用格式为： Dy=diff(Y)：计算向量Y的向前差分，并把结果赋值给向量Dy Dy(i)=Y(i+1)-Y(i)，i=1,2,…,n-1。注意向量Dy元素个数比Y少 Dy=diff(Y,n)：计算向量Y的n阶向前差分。例如， diff(Y,2)=diff(diff(Y))=DX(i+1)-DX(i)= Y(i+2)-2Y(i+1)+Y(i) ， i=1,2 ……n-2。 DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。 A=pascal(4) A = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 B=diff(A) B = 0 1 2 3 0 1 3 6 0 1 4 10 C=diff(A,2) C = 0 0 1 3 0 0 1