第3章 matlab微积分问题的计算机求解-医学课件.ppt

建立用拟合（插值）多项式计算各阶导数的poly_drv.m function der=poly_drv(xd,yd,a) m=length(xd)-1; d=polyfit(xd-a,yd,m); c=d(m:-1:1); ％去掉常数项 fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);end der=c.*fact; 例： xd=[ 0 0.2000 0.4000 0.6000 0.8000 1.000]; yd=[0.3927 0.5672 0.6982 0.7941 0.8614 0.9053]; a=0.3; der=poly_drv(xd,yd,a) der = 0.6533 -0.9710 1.0406 -1.3750 1.8750 * 3.3.4 二元函数的梯度计算 格式： 若z矩阵是建立在等间距的形式生成的网格基础上，则实际梯度为 * 例： 计算梯度，绘制引力线图： [x,y]=meshgrid(-3:.2:3,-2:.2:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); [fx,fy]=gradient(z); fx=fx/0.2; fy=fy/0.2; contour(x,y,z,30); hold on; quiver(x,y,fx,fy) %绘制等高线与 引力线图 * 绘制误差曲面: zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y); zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y); surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.08]) figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.11]) ％建立一个新图形窗口 * 为减少误差，对网格加密一倍： [x,y]=meshgrid(-3:.1:3,-2:.1:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); [fx,fy]=gradient(z); fx=fx/0.1; fy=fy/0.1; zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y); zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y); surf(x,y,abs(fx-zx)); axis([-3 3 -2 2 0,0.02]) figure; surf(x,y,abs(fy-zy)); axis([-3 3 -2 2 0,0.06]) * 3.4 数值积分问题 4.3.1 由给定数据进行梯形求积 Sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2 * 格式： S=trapz(x,y) 例： x1=[0:pi/30:pi]; y=[sin(x1) cos(x1) sin(x1/2)]; x=[x1 x1 x1]; S=sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2 S = 1.9982 0.0000 1.9995 S1=trapz(x1,y) % 得出和上述完全一致的结果 S1 = 1.9982 0.0000 1.9995 * 例： 画图 x=[0:0.01:3*pi/2, 3*pi/2]; % 这样赋值能确保 3*pi/2点被包含在内 y=cos(15*x); plot(x,y) % 求取理论值 syms x, A=int(cos… (15*x),0,3*pi/2) A = 1/15 * 随着步距h的减小，计算精度逐渐增加： h0=[0.1,0.01,0.001,0.0001,0.00001,0.000001]; v=[]; for h=h0, x=[0:h:3*pi/2, 3*pi/2]; y=cos(15*x); I=trapz(x,y); v=[v; h, I, 1/15-I ]; end v v = 0.1000 0.0539 0.0128 0.0100 0.0665 0.0001 0.0010 0.0667 0.0000 0.0001 0.0667 0.0000 0.0000 0.0667 0.0000