matlab数值分析程序--高等数学,数值代数的matlab实现-文字版, matlab电子书, 和matlab 有关的电子书:

9.8.2 二重积分的梯形公式及其MATLAB程序

9.8.2  二重积分的梯形公式及其MATLAB程序

(一) 二重积分的基本梯形公式及其MATLAB程序

例9.8.3  根据基本梯形公式编写MATLAB程序,计算d近似值,并将计算结果与精确值比较.其中Dxy是矩形区域.

解 编写并输入下列MATLAB程序

>> a=0;b=4;c=-1;d=2; fac=sin(a+c)/(a+c);fad=sin(a+d)/(a+d);

fbc=sin(b+c)/(b+c);fbd=sin(b+d)/(b+d); h=(b-a)*(d-c);

F=fac+fad+fbc+fbd; T2=(h*F)/4

syms x y

bjh=sin(x+y)/(x+y); jfx=int(bjh,x,a,b);

jfy=int(jfx,y,c,d); I2=double(jfy), juewu=abs(I2-T2)

运行后屏幕显示如下

T2 =                    I2 =                 juewu =

  3.88977135362262        3.64403676611252         0.24573458751009

 

 (二) 二重积分的复合梯形公式及其MATLAB程序

例9.8.4  根据(9.93)式编写MATLAB程序计算d近似值,分别取,并将计算结果与精确值比较.观察节点的个数与误差的关系.问取时,计算机能计算吗,为什么?其中Dxy是矩形区域.

解 (1)编写下列MATLAB主程序并保存名为r2trapf2.m的M文件

function [T2,I2,juewu]=r2trapf2(a,b,n,c,d,m)

x=zeros(1,n);y=zeros(1,m); P=zeros(n,m); fyx=zeros(n,m);

hx=(b-a)/(n-1);hy=(d-c)/(m-1); x(1)=a;x(n)=b;y(1)=c;y(m)=d;

for k=1:n-1

x(k+1)=a+k*hx;

for j=1:m-1

y(j+1)=c+j*hy;

fyx(k+1,j+1)=exp(-(x(k+1)^2+y(j)^2))/(x(k+1)+y(j));

fyx(1,j+1)=exp(-(x(1)^2+y(j+1)^2))/(x(1)+y(j+1));

fyx(k+1,1)=exp(-(x(k+1)^2+y(1)^2))/(x(k+1)+y(1));

fyx(1,1)=exp(-(x(1)^2+y(1)^2))/(x(1)+y(1));

 fyx(1,m)=exp(-(x(1)^2+y(m)^2))/(x(1)+y(m));

fyx(n,1)=exp(-(x(n)^2+y(1)^2))/(x(n)+y(1));

fyx(n,m)=exp(-(x(n)^2+y(m)^2))/(x(n)+y(m));

P(k+1,j+1)=4;P(1,j+1)=2;P(n,j+1)=2;

P(k+1,1)=2; P(k+1,m)=2;

P(1,1)=1; P(1,m)=1; P(n,1)=1; P(n,m)=1;

end

end

P;fyx; F=P.*fyx; F1=sum(F); F2=sum(F1); T2=hx*hy*F2/4;

syms t w

bjh=exp(-(w^2+t^2))/(w+t); jfx=int(bjh,w,a,b);

jfy=int(jfx,t,c,d); I2=double(jfy); juewu=abs(I2-T2);

(2)输入下列MATLAB程序

>>a=0;b=4;n2=2;c=1;d=2;

m2=2;[T2,I2,juewu2]=r2trapf2(a,b,n2,c,d,m2)

n12=12;;m12=12;

[T12,I12,juewu12]=r2trapf2(a,b,n2,c,d,m12)

n112=112;;m112=112;

[T112,I2,juewu112]=r2trapf2(a,b,n112,c,d,m112)

n1112=1112;

m1112=1112;

[T1112,I2,juewu1112]=r2trapf2(a,b,n1112,c,d,m1112)

运行后屏幕显示如下

Warning: Explicit integral could not be found.

T2 =                  I2 =                  juewu2 =

   0.37703726923921     0.06889208904324          0.30814518019597

T12 =                juewu12 =              T112 =

  0.21706671440746      0.14817462536422          0.07079360039872

,juewu112 =           T1112 =               ,juewu1112 =

0.00190151135548      0.06908616401326      1.940749700208577e-004

由此可见,节点的个数越多,近似值与精确值的绝对误差越小.但是,取时,计算机就不能计算了,因为计算机内存不够.见下例:

>> a=0;b=4;n2=111112;c=1;d=2;m2=111112;

[T2,I2,juewu2]=r2trapf2(a,b,n2,c,d,m2)

??? Error using ==> zeros

Product of dimensions is greater than maximum integer.

 

Error in ==> D:\MATLAB6p5\work\r2trapf2.m

On line 2  ==> x=zeros(1,n);y=zeros(1,m);

P=zeros(n,m); fyx=zeros(n,m);

用复合梯形公式(9.93)数值计算二重积分I2=d的MATLAB主程序

function [T2,I2,juewu]=r2trapf2(fun,a,b,n,c,d,m)

x=zeros(1,n);y=zeros(1,m); P=zeros(n,m); fyx=zeros(n,m);

hx=(b-a)/(n-1);hy=(d-c)/(m-1); x(1)=a;x(n)=b;y(1)=c;y(m)=d;

for k=1:n-1

x(k+1)=a+k*hx;

for j=1:m

y(j+1)=c+j*hy; fyx(k+1,j+1)=feval(fun,x(k+1),y(j));

fyx(1,j+1)=feval(fun,x(1),y(j+1));

fyx(k+1,1)=feval(fun,x(k+1),y(1));

fyx(1,1)=feval(fun,x(1),y(1));

 fyx(1,m)=feval(fun,x(1),y(m));

fyx(n,1)=feval(fun,x(n),y(1));

 fyx(n,m)=feval(fun,x(n),y(m));

P(k+1,j+1)=4; P(1,j+1)=2;P(n,j+1)=2; P(k+1,1)=2;

P(k+1,m)=2;

P(1,1)=1; P(1,m)=1; P(n,1)=1; P(n,m)=1;

end

end

P;fyx; F=P.*fyx; F1=sum(F); F2=sum(F1); T2=hx*hy*F2/4;

syms t w

bjh=feval(fun,w,t); jfx=int(bjh,w,a,b);

 jfy=int(jfx,t,c,d);

I2=double(jfy); juewu=abs(I2-T2);

欢迎转载,转载请注明来自一手册:http://yishouce.com/book/3/3090802.html
友情链接It题库(ittiku.com)| 版权归yishouce.com所有| 友链等可联系 admin#yishouce.com|粤ICP备16001685号-1