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

9.7.2 无穷积分的近似计算及其MATLAB程序

9.7.2  无穷积分的近似计算及其MATLAB程序

 

 (一) 累积求和无穷积分法及其MATLAB程序

累积求和无穷积分法计算无穷积分的近似值的MATLAB主程序

function [k,suj,wugu,Fjj,WC]=wqjfjx(a,Wucha,m)

r=a;k=0; wugu=1;

syms t

Y=input( '计算没有结束,请输入被积函数f(x)=');%exp(-t)./(1+t.^4);

su=int(Y,t,a,a+2)

suj= double(su)

while ((wugu>Wucha)&(k<m))

k=k+1,r=a+2^k,r1=a+2^(k+1)

syms t

Y=input( '计算没有结束,请输入被积函数f(x)=');

F1=int(Y,t,r,r1), intF1= double(F1);

suj=suj+intF1,wugu=abs(intF1)

Fj=int(Y,t,0,+inf)

Fjj= double(Fj), WC=abs(Fjj- suj)

end

disp('k是累加次数,suj是利用累积求和无穷积分法计算f(x)从a到正无穷的无穷积分的近似值,wugu是f(x)在[r(n),r(n+1)]上的积分的绝对值,Fjj是f(x)从a到正无穷的无穷积分的值,WC=|Fjj-suj|')

 

例9.7.4   计算d,使精确度为,并与精确值比较.

解 取,根据(9.72)式计算d.输入下面的程序:

>> a=0; Wucha =1.e-6;m=100;

[k,suj,wugu,Fjj,WC]=wqjfjx(a,Wucha,m)

运行后屏幕显示

计算没有结束,请输入被积函数f(x)=

输入函数:exp(-t)./(1+t.^4);

运行后屏幕显示

Warning: Explicit integral could not be found.

> In C:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58

  In C:\MATLAB6p5\work\wqjfjx.m at line 5

su =

int(exp(-t)/(1+t^4),t = 0 .. 2)

suj =

   0.62745967782312

k =

     1

r =

     2

r1 =

     4

计算没有结束,请输入被积函数f(x)=

输入函数:exp(-t)./(1+t.^4);

运行后屏幕显示提示,……,如此继续下去,当wugu<时,屏幕显示如下内容:

k是累加次数,suj是利用累积求和无穷积分法计算f(x)从a到正无穷的无穷积分的近似值,wugu是f(x)在[r(n),r(n+1)]上的积分的绝对值,Fjj是f(x)从a到正无穷的无穷积分的值,WC=|Fjj-suj|

k =

     3

suj =

   0.63047783491711

wugu =

    5.598881094648278e-008

Fjj =

   0.63047783491850

WC =

    1.386668557756821e-012

 

(二) 无穷区间的截断法

 

例9.7.5   用数值方法计算ed,使精确度为,并与精确值比较.

解  为估计截断部分ed的值,利用

eded.

将下面的函数保存为名为nm.m的M文件

function y=nm(x)

y=exp(-3*x.^2)./(3*x);

输入程序

>> N=2:5;y=nm(N)

运行后屏幕显示如下

y =  1.0e-005 *

0.10240353922214  0.00000002088365  0.00000000000000  0.00000000000000

由此可见,如果不计算误差传播,取,可满足ed.

输入程序

>> x=0:0.1:3; y=exp(-3*x.^2); z=trapz(x,y)

syms t,

F1=int(exp(-3*t.^2),t,0,+inf),F3=int(exp(-3*t.^2),t,0,3)

intF1= double(F1), Wuchaz=abs(intF1-z)

intF3= double(F3), Wucha3=abs(intF1-intF3)

运行后屏幕显示当时,分别用梯形公式和int函数计算ed的近似值z和精确解 F3(或intF3)及其与的误差Wuchaz和Wucha3如下

z =                      F1 =                    intF1 = 0.51166335397311        1/6*3^(1/2)*pi^(1/2)      0.51166335397324

F3 =                                             intF3 =

1/6*erf(3*3^(1/2))*3^(1/2)*pi^(1/2)       0.51166335397314

Wuchaz =                                 Wucha3 =

    1.294520046712933e-013               1.025846074753645e-013

 

例9.7.6   如果取,根据(9.72)式,当取何值时, d的近似值的截断误差为,计算其近似值,并与精确值比较.

解 (1)首先求的值.输入下面的程序

>> k=ceil(sqrt(10^6/(4*pi)))

运行后屏幕显示

k =

283

因此,当时,的近似值的截断误差为.

(2)然后计算的近似值,使截断误差为,并与精确值比较.

编写并保存名为wqjfjx1.m的计算主程序如下

function [k,suj,wugu]=wqjfjx1(a,m)

r=a;k=0; wugu=1; syms t

fun=sin(t)./(1+t.^2); su=int(fun,t,a,2*pi);

suj= double(su);

while (k<m)

k=k+1,r=2*k*pi; r1=2*(k+1)*pi;

syms t

F1=int(fun,t,r,r1); intF1= double(F1);

suj=suj+intF1, wugu=abs(intF1),

end

输入下面的程序

>>  a=0,m=283,[k,suj,wugu]=wqjfjx1(a,m)

syms t

fun= sin(t)./(1+t.^2); Fj=int(fun,t,0,+inf),

Fjj= double(Fj)

WC=abs(Fjj- suj)

运行后屏幕显示计算过程和结果如下:

k =   suj =                                    wugu =

  1    0.64068026831490 - 0.00000000000000i       0.01602703309687

k =   suj =                                    wugu =

  2    0.64399926393154 - 0.00000000000000i       0.00331899561664

k =   suj =                                    wugu =

  3    0.64519499777913 - 0.00000000000000i       0.00119573384758

…………………………………………………………

k =    suj =                                   wugu =

283  0.64676080872636 - 0.00000000000000 i  2.223370580497967e-009

Fj =

-i*sinint(i)*cosh(1)-cosint(i)*sinh(1)+1/2*i*sinh(1)*pi

Fjj =                                         WC =

  0.64676112277913 - 0.00000000000000i      3.140527699008189e-007

 

(三) 高斯-拉盖尔求积公式及其MATLAB程序

计算次正交多项式——拉盖尔多项式和系数公式的MATLAB主程序

function [K,UK,AK]=laguerre(n)

syms x

K=0,U0=1,K=1,U1=1-x,u1=1;u2=1-x;

duk=diff(u2,x); n=n-1;sun1=1;

duk2= duk^2;

A1=(sun1^2)/(x*duk2);A1=simple(A1),

sun=1;

for k=1:2:n

   u3=(1+2*k-x)*u2-(k^2)*u1;K=k+1, UK =simple(u3),

   duk1=diff(u3,x);

   sun2=sun1*K;

duk21= duk1^2;

   Ak1=(sun1^2)/(x*duk21); AK =simple(Ak1),

   u4=(1+2*(k+1)-x)*u3-((k+1)^2)*u2;K=k+2, UK =simple(u4),

   duk=diff(u4,x);

sun3=sun2*K;sun1=sun3;

duk2= duk^2;

    Ak2=(sun3^2)/(x*duk2); AK =simple(Ak2),

    u1=u3;u2=u4;

if k>n,return,end;

end

K=k+2; UK; AK;

 

例9.7.7  根据(9.73)式、(9.74)式和(9.75)式,计算1至7次正交多项

式——拉盖尔多项式、节点的横坐标和系数公式及其在处的值.

解 (1)建立并保存以laguerre.m文件命名的M文件函数.

(2)计算1至7次正交多项式——拉盖尔多项式和系数公式.输入程序

>>n=7,[K,U,Ak]=laguerre(n)

运行后输出结果

K =

     0

U0 =

     1

K =

     1

U1 =

1-x

A1 =

1/x

K =

     2

UK =

2-4*x+x^2

AK =

1/x/(-4+2*x)^2

K =

     3

UK =

6-18*x+9*x^2-x^3

AK =

4/x/(6-6*x+x^2)^2

K =

     4

UK =

-16*x^3+72*x^2-96*x+24+x^4

AK =

9/4/x/(-24+36*x-12*x^2+x^3)^2

K =

     5

UK =

-200*x^3+600*x^2-600*x+120+25*x^4-x^5

AK =

576/x/(120*x^2-240*x+120-20*x^3+x^4)^2

K =

     6

UK =

720-4320*x+5400*x^2-2400*x^3+450*x^4-36*x^5+x^6

AK =

400/x/(-720+1800*x-1200*x^2+300*x^3-30*x^4+x^5)^2

K =

     7

UK =

5040-35280*x+52920*x^2-29400*x^3+7350*x^4-882*x^5+49*x^6-x^7

AK =

518400/x/(-15120*x+5040-4200*x^3+12600*x^2-42*x^5+630*x^4+x^6)^2

(3)求1至7次正交多项式——拉盖尔多项式的零点,即求积节点. 输入程序

>>syms x y

y1=('1-x=0')

y2=('x^2-4*x+2=0')

y3=('-x^3+9*x^2-18*x+6=0')

y4=('x^4-16*x^3+72*x^2-96*x+24=0')

y5=('-x^5+25*x^4-200*x^3+600*x^2-600*x+120=0')

y6=('720-4320*x+5400*x^2-2400*x^3+450*x^4-36*x^5+x^6=0')

y7=('5040-35280*x+52920*x^2-29400*x^3+7350*x^4-882*x^5+49*x^6-x^7=0')

y11=solve(y1,x), x1=double(y11)

y22=solve(y2,x) , x2=double(y22),

y33=solve(y3,x) , x3=double(y33),

y44=solve(y4,x) , x4=double(y44),

y55=solve(y5,x) , x5=double(y55),

y66=solve(y6,x) , x6=double(y66),

y77=solve(y7,x) , x7=double(y77),

运行后可以得到(9.74)式中的节点的横坐标(略).

(4)求在节点的横坐标处的系数. 输入程序

>> x1=[1],A1 =1./x1

x2=[3.41421356237309,0.58578643762690],A2 =1./x2./(-4+2*x2).^2

x3=[6.28994508293748 0.41577455678348 2.29428036027904],

 A3=4./x3./(6-6*x3+x3.^2).^2

x4=[9.39507091230113 4.53662029692113 1.74576110115835 0.32254768961939],

A4 =9/4./x4./(-24+36*x4-12*x4.^2+x4.^3).^2

x5=[0.26356031971814 1.41340305910652 3.59642577104072 7.08581000585884 12.64080084427578],

A5 =576./x5./(120-240*x5+120*x5.^2-20*x5.^3+x5.^4).^2

x6=[0.22284660417926 1.18893210167262 2.99273632605931 5.77514356910451 9.83746741838259 15.98287398060170],

 A6 =400./x6./(-720+1800*x6-1200*x6.^2+300*x6.^3-30*x6.^4+x6.^5).^2

x7=[0.19304367656036 1.02666489533919 2.56787674495075 4.90035308452648 8.18215344456286 12.73418029179782 19.39572786226254],

A7 =518400./x7./(5040-15120*x7+12600*x7.^2-4200*x7.^3+630*x7.^4-42*x7.^5+x7.^6).^2

运行后输出节点的横坐标和对应的系数AK (K如下

x1 =

     1

A1 =

     1

x2 =

   3.41421356237309   0.58578643762690

A2 =

   0.03661165235168   0.21338834764832

x3 =

   6.28994508293748   0.41577455678348   2.29428036027904

A3 =

   0.01038925650159   0.71109300992917   0.27851773356924

x4 =

  9.39507091230113   4.53662029692113   1.74576110115835   0.32254768961939

A4 =

 0.00003370591910   0.00243049428219   0.02233866827736   0.03769713152135

x5 =

  Columns 1 through 4

 0.26356031971814   1.41340305910652   3.59642577104072   7.08581000585884

  Column 5

 12.64080084427578

A5 =

  Columns 1 through 4

 0.52175561058281   0.39866681108317   0.07594244968171   0.00361175867992

  Column 5

 0.00002336997239

x6 =

  Columns 1 through 4

 0.22284660417926   1.18893210167262   2.99273632605931   5.77514356910451

  Columns 5 through 6

 9.83746741838259  15.98287398060170

A6 =

  Columns 1 through 4

 0.01274901872083   0.01158335641034   0.00314926061317   0.00028886659592

  Columns 5 through 6

 0.00000725047786   0.00000002495966

x7 =

  Columns 1 through 4

 0.19304367656036   1.02666489533919   2.56787674495075   4.90035308452648

  Columns 5 through 7

 8.18215344456286  12.73418029179782  19.39572786226254

A7 =

  Columns 1 through 4

   0.40931895170127   0.42183127786172   0.14712634865750   0.02063351446872

  Columns 5 through 7

   0.00107401014328   0.00001586546435   0.00000003170315

用高斯—拉盖尔无穷积分公式计算ed的数值积分及其截断误差公式的MATLAB主程序

function [GL,Y, RGn]=GaussL1(fun,X,A)

n=length(X);n2=2*n; Y=feval(fun,X);

GL=sum(A.*Y); sun=1; su2n=1; su2n1=1;

for k=1:n

sun=sun*k;

end

for k=1:n2

su2n=su2n*k;

end

syms M

RGn= ( sun^2)*M/( su2n);

用高斯—拉盖尔无穷积分公式计算ed的数值积分和误差估计的MATLAB主程序

function [GL,Y,Rn]=GaussL2 (fun,X,A,fun2n)

n=length(X);n2=2*n; Y=feval(fun,X);

GL=sum(A.*Y); sun=1; su2n=1; su2n1=1;

for k=1:n

sun=sun*k;

end

for k=1:n2

su2n=su2n*k;

end

mfun2n1 =max(fun2n);

mfun2n=abs(mfun2n1);

Rn= ( sun^2)*mfun2n /( su2n);

 

例9.7.8  用高斯—拉盖尔无穷积分公式计算d,取n=7,再根据截断误差公式写出误差公式,并将计算结果与精确值进行比较.

解 (1)建立并保存以fun.m文件命名的M文件函数

function y = fun(x)

       y = 1./(5+x.^2);

(2)建立并保存用GaussL2.m文件命名的M文件函数.

(3)输入程序

>> syms x

fun2n=diff(exp(-x)./(5+x^2),x,10); fun2ns=simple(fun2n)

    (4)运行后,再输入程序

>>X=[0.19304367656036,1.02666489533919,2.56787674495075, 4.90035308452648, 8.18215344456286,12.73418029179782,19.39572786226254];

A=[0.40931895170127,0.42183127786172,0.14712634865750, 0.02063351446872,0.00107401014328,0.00001586546435,0.00000003170315];

x=0:0.001:1000;

fun2ns=exp(-x).*(50889062500*x-29608750000*x.^2+89352450000*x.^5+69284490000*x.^6+9796290000*x.^7-185554687500*x.^3-70161046875*x.^4-6677396250*x.^8-868165200*x.^10-3734325000*x.^9+14759850*x.^12+1748400*x.^14-79569000*x.^11+7568400*x.^13+36675*x.^16+320*x.^18+x.^20+285840*x.^15+3780*x.^17+20*x.^19+7000234375)./(5+x.^2).^11;

[GL,Y, Rn]=GaussL2(@fun,X,A,fun2ns)

syms t

fi=int(exp(-t)/(5+t^2),t,0,+inf); Fs= double (fi), wGL = double (abs(fi-GL))

运行后屏幕显示如下

GL =

   0.16435187219536

Y =

  Columns 1 through 4

0.19852039332099  0.16517893285268  0.08625157652455  0.03446676087003

  Columns 5 through 7

 0.01389899751582   0.00598231575880   0.00262333436453

Rn =

   0.06539672874780

Fs =

  0.16435041671009 + 0.00000000000000i

wGL =

    1.455485269162850e-006

 

(四)高斯—埃尔米特求积公式及其MATLAB程序

高斯—埃尔米特无穷数值积分公式计算ed及其截断误差公式的MATLAB主程序

function [GH,Y,RHn]=GaussH1(fun,X,A)

n=length(X);n2=2*n; Y=feval(fun,X); GH=sum(A.*Y);

sun=1; su2n=1; su2n1=1;

for k=1:n

sun=sun*k;

end

for k=1:n2

su2n=su2n*k;

end

syms M

RHn = ( sun*sqrt(pi))*M/( 2^n*su2n);

例9.7.9  用高斯—埃尔米特无穷积分公式计算ed,取n=5,再根据截断误差公式写出误差公式,并将计算结果与精确值进行比较.

解 (1)建立并保存fun.m文件命名的M文件函数

function y = fun (x)

              y = sin(3*x)./sqrt(2*pi);

(2)建立并保存GaussH1.m文件命名的M文件函数.

(3)根据表9–8 ,输入程序

>> X=[2.0201828705 -2.0201828705 0.9585724646 -0.9585724646 0];

A=[0.0199532421  0.0199532421  0.3936193232  0.3936193232  0.9453087205];

[GH,Y, RHn]=GaussH1(@fun,X,A)

syms x

fi1=int(exp(-x^2).*sin(3*x)./(sqrt(2*pi)),x,0,+inf);

fi2=int(exp(-x^2).*sin(3*x)./(sqrt(2*pi)),x, -inf,0);

fi=fi1+fi2; Fs= double (fi), wGL = double (abs(fi-GH))

运行后屏幕显示取n=5,用高斯—埃尔米特无穷积分公式数值计算的结果GH、截断误差公式RHn和精确值Fs及其两者的绝对误差wGL依次如下:

GH =

     0

Y =

  Columns 1 through 4

-0.08808725583239  0.08808725583239  0.10482362982925  -0.10482362982925

  Column 5

                  0

RHn =

2494507032021713/1361888527316837990400*M

Fs =

     0

wGL =

     0

高斯—埃尔米特无穷数值积分公式计算ed及其误差估计的MATLAB主程序

function [GH,Y, RHn]=GaussH2(fun,X,A,fun2n)

n=length(X);n2=2*n; Y=feval(fun,X);

GH=sum(A.*Y); sun=1; su2n=1; su2n1=1;

for k=1:n

sun=sun*k;

end

for k=1:n2

su2n=su2n*k;

end

mfun2n =max(fun2n);

RHn = ( sun*sqrt(pi))* mfun2n /( 2^n*su2n);

 

例9.7.10  用高斯—埃尔米特无穷积分公式计算d,取n=5,再根据截断误差公式写出误差公式,并将计算结果与精确值进行比较.

解 (1)建立并保存fun.m文件命名的M文件函数

function y = fun(x)

       y = 1./(5+x.^2);

(2)建立并保存GaussH2.m文件命名的M文件函数.

(3)输入程序

>> syms x

fun2n=diff(exp(-x.^2)./(5+x.^2),x,10); fun2ns=simple(fun2n)

    (4)运行后再输入程序

>> X=[2.0201828705 -2.0201828705 0.9585724646 -0.9585724646 0];

A=[0.0199532421 0.0199532421 0.3936193232 0.3936193232 0.9453087205];

x=-1000:0.001:1000;

fun2ns=32*exp(-x.^2).*(425471484375*x.^2+26348013150*x.^10-35644522500*x.^8-290495953125*x.^4-262342473750*x.^6+13793263150*x.^12+2639594700*x.^14+34752150*x.^16-88200525*x.^18-18517905*x.^20-1466910*x.^22+32520*x.^24+15840*x.^26+1200*x.^28+32*x.^30-32395781250)./(5+x.^2).^11;

[GH,Y, RHn]=GaussH2(@fun,X,A, fun2ns)

syms t

fi1=int(exp(-t.^2)./(5+t.^2),t,0,+inf); fi2=int(exp(-t.^2)./(5+t.^2),t, -inf,0);

fi=fi1+fi2; Fs= double (fi), wGL = double (abs(fi-GH))

运行后屏幕显示如下

GH =

   0.32646126745168

Y =

  Columns 1 through 4

0.11011834734512 0.11011834734512 0.16895142009535 0.16895142009535

  Column 5

   0.20000000000000

RHn =

   0.02897468774679

Fs =

   0.32640983502878

wGL =

    5.143242289669318e-005

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