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

9.6.4 拉道积分公式和洛巴托积分公式及其MATLAB程序

9.6.4  拉道积分公式和洛巴托积分公式及其MATLAB程序

 

    (一)  用洛巴托数值积分的MATLAB程序

调用格式一:Q = quadl (FUN,a,b)

调用格式二:Q = quadl (FUN,a,b,TOL)

调用格式三:[Q,FCNT] = quadl (...)

调用格式四:quadl (FUN,a,b,TOL,TRACE)

调用格式四:quadl (FUN,A,B,TOL,TRACE,P1,P2,...)

调用格式五:quadl FUN(X,P1,P2,...)

例9.6.7  用函数quadl和quad8及quad分别数值计算ed,分别取精度为,将计算结果分别与精确值比较,并画出图形.

解  输入程序

>>  subplot(2,1,1),

[Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-6)

[Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-6,3)

[Q,FCNT]=quad(@fun,0,pi/2,1.e-6)

syms x

fi=int(exp((-x.^2)./2)./(sqrt(2*pi)),x,0, pi/2);

Fs=double(fi),

wQl=double(abs(fi-Ql))

wQ8=double(abs(fi-Q8)),

wQ=double(abs(fi-Q))

xlabel('自变量 X'), ylabel('因变量 Y')

title('取精度为1.e-6,被积函数f(x)=exp(-x^2/2)/(sqrt(2 pi))在[0,pi/2]上的数值积分')

legend('取精度1.e-6,用牛顿–科茨公式计算积分过程的散点图')

hold on

subplot(2,1,2)

[Ql14,FCNTl14]=quadl(@fun,0,pi/2,1.e-14)

[Q814,FCNT814]=quad8(@fun,0,pi/2,1.e-14,3)

[Q14,FCNT14]=quad(@fun,0,pi/2,1.e-14,1)

wQl14=double(abs(fi-Ql14)),

wQ814=double(abs(fi-Q814))

wQ14=double(abs(fi-Q14)),

xlabel('自变量 X'), ylabel('因变量 Y')

title('取精度为1.e-14,被积函数f(x)=exp(-x^2/2)/(sqrt(2 pi))在[0,pi/2]上的数值积分')

legend('取精度1.e-14,用牛顿–科茨公式计算积分过程的散点图')

hold off

运行后屏幕显示如下

Ql =                   FCNTl =     Q8 =                   FCNT8 =

0.44188501970377         18        0.44188501721659        33

Q =                    FCNT =       Fs =

0.44188511819589         13        0.44188501721659

wQl =                               wQ8 =

2.487183090649321e-009          3.572753311048285e-015

wQ =                               Ql14 =

1.009793004969618e-007          0.44188501721659,

FCNTl14 =       Q814 =                      FCNT814 =

138          0.44188501721659              129

Q14 =                   FCNT14 =       wQl14 =

0.44188501721659          553        1.310619347103004e-016

wQ814 =                            wQ14 =

9.098267021473094e-017             2.003963224778471e-017

 

例9.6.8  用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算ed,取精度分别为,并与精确值比较.请再改用取10个节点,分别用三阶拉格朗日插值多项式和三次样条构造被积函数,然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算,取计算精度分别为.再改用100, 1 000,10 000个节点,用三次样条构造被积函数计算,观察计算结果与精确值的绝对误差有何变化.

解  (1)建立名为fun.m的M文件

function y = fun(x)

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

取精度为,输入程序

>> F=inline('exp(x.^3).*cos(2*x)'); [Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-4,3)

[Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-4,3)

[Q,FCNT]=quad(@fun,0,pi/2,1.e-4,3)

syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2);

Fs=double(fi), wQl=double(abs(fi-Ql) )

wQ8=double(abs(fi-Q8) ), wQ=double(abs(fi-Q) )

运行后屏幕显示如下

Ql =

  -6.45505792320975

FCNTl =

    48

Q8 =

  -6.45505883179317

FCNT8 =

    33

Q =

  -6.45506213501363

FCNT =

    33

Warning: Explicit integral could not be found.

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

Fs =

  -6.45505792116786

wQl =

    2.041896319084332e-009

wQ8 =

    9.106253137107240e-007

wQ =

    4.213845774426019e-006

(2)取精度为,输入程序

>> F=inline('exp(x.^3).*cos(2*x)');

[Ql,FCNTl]=quadl(@fun,0,pi/2,1.e-14,3)

[Q8,FCNT8]=quad8(@fun,0,pi/2,1.e-14,3)

[Q,FCNT]=quad(@fun,0,pi/2,1.e-14,3)

syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2);

Fs=double(fi), wQl=double(abs(fi-Ql))

wQ8=double(abs(fi-Q8)), wQ=double(abs(fi-Q))

运行后屏幕显示如下

Ql =

  -6.45505792116785

FCNTl =

   798

Q8 =

  -6.45505792116785

FCNT8 =

        4273

Q =

  -6.45505792116785

FCNT =

        2961

Warning: Explicit integral could not be found.

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

Fs =

  -6.45505792116786

wQl =

    3.663345763769995e-015

wQ8 =

    1.886988924369744e-015

wQ =

    2.775167344069869e-015

这三种方法计算结果的绝对误差都小于,达到所要求的精度.

(3)取10个节点,改用三阶拉格朗日插值多项式构造被积函数,然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算,取精度分别为,输入程序

>>X=[0,0.1745,0.3491,0.5236,0.6981,0.8727,1.0472,1.2217,1.3963,1.5708];

Y= exp(X.^3).*cos(2* X); L1=lfun(X(1:4),Y(1:4)),

L2=lfun(X(4:7),Y(4:7)), L3=lfun(X(7:10),Y(7:10)),

运行后输出多项式L1,L2,L3.然后输入程序

 >> syms x

L1 = inline(' 528475254238581253770784020141/1267650600228229401496703205376*x.^3-2150921733462980567549357641495/1267650600228229401496703205376*x.^2-110805035875068143559156492912389/3318709271397504573118368991674368*x+1');

L2 = inline(' -32514595309399085335529842115599/2535301200456458802993406410752*x.^3+28880777683589964816385504370251/1267650600228229401496703205376*x.^2-9690424145215092365005571145155/633825300114114700748351602688*x+21181277275320971857213583641339/5070602400912917605986812821504');

L3 = inline(' -22347657567307906685985592460479/39614081257132168796771975168*x.^3+19439764627825105995989937383455/9903520314283042199192993792*x.^2-45346855519117547561405281365707/19807040628566084398385987584*x+17651524197374204698139309247929/19807040628566084398385987584');

[Q1,FCNT1] = quad(L1,0, 0.5236,1.e-4);

[Q2,FCNT2] =quad(L2, 0.5236, 1.0472,1.e-4);

[Q3,FCNT3] = quad(L3, 1.0472, pi/2,1.e-4);

Q= Q1+ Q2+ Q3 ,FCNT= FCNT1+ FCNT2+ FCNT3

[Q81,FCNT81] = quad8(L1,0, 0.5236,1.e-4);

[Q82,FCNT82] =quad8(L2, 0.5236, 1.0472,1.e-4);

[Q83,FCNT83] = quad8(L3, 1.0472, pi/2,1.e-4);

Q8= Q81+ Q82+ Q83 ,FCN8T= FCNT81+ FCNT82+ FCNT83

[Ql1,FCNTl1] = quadl(L1,0, 0.5236,1.e-4);

[Ql2,FCNTl2] =quadl(L2, 0.5236, 1.0472,1.e-4);

[Ql3,FCNTl3] = quadl(L3, 1.0472, pi/2,1.e-4);

Ql= Ql1+ Ql2+ Ql3 ,FCNTl= FCNTl1+ FCNTl2+ FCNTl3

syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2);

Fs= double (fi), wQl= double (abs(fi-Ql) )

wQ8= double (abs(fi-Q8) ), wQ= double (abs(fi-Q) )

运行后屏幕显示如下

Q =                     FCNT=    Q8 =                FCN8T =

-6.63556515894584        39      -6.63556515894593         99

Ql =                     FCNTl =

-6.63556515894607            54

Warning: Explicit integral could not be found.

> In D:\MATLAB6P5\toolbox\symbolic\@sym\int.m at line 58

Fs =                            wQl =

-6.45505792116786             0.18050723777821

wQ8 =                           wQ =

0.18050723777807              0.18050723777798

(4)取10个节点,改用三次样条构造被积函数,然后用辛普森求积公式、牛顿–科茨求积公式和洛巴托求积公式分别计算,取精度分别为.

首先建立并保存名为scjfQl8.m的M文件

function wQ=scjfQl8(n)

a=0; b= pi/2;X=0:pi/(2*n-1):pi/2; fs= exp(X.^3).*cos(2* X);

pp=spline(X,fs); Qsl=quadl(@ppval,a,b,[],[],pp)

Qs8=quad8(@ppval,a,b,[],[],pp),

Qs=quad(@ppval,a,b,[],[],pp)

syms x

fi=int(exp(x.^3).*cos(2*x),x,0, pi/2); Fs=double(fi)

wQl=double(abs(fi-Qsl)),

wQ8=double(abs(fi-Qs8))

wQ=double(abs(fi-Qs))

然后输入程序

>> n=10;wQ=scjfQl8(n)

运行后屏幕显示结果如下

Qsl =                 Qs8 =                   Qs =

  -6.33658113206801     -6.33658632890273    -6.33659053214403

Warning: Explicit integral could not be found.

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

  In C:\MATLAB6p5\work\scjfQl8.m at line 7

Fs =

  -6.45505792116786

wQl =                 wQ8 =                   wQ =

0.11847678909985      0.11847159226513      0.11846738902383

 

例9.6.10  测得定积分d中被积函数在积分变量的某几个特定X =0, 0.157 1, 0.471 2, 0.549 8, 0.628 3, 0.863 9, 1.021 0, 1.099 6, 1.570 8处的值分别为Y= 0, 0.156 5, 0.454 0, 0.522 5, 0.587 8, 0.760 4, 0.852 6, 0.891 0, 1.000 0,试根据这些值首先分别用分段二阶、三阶拉格朗日插值多项式和三次样条插值构造被积函数,然后分别用辛普森和洛巴托求积公式计算,取精度分别为,并与精确值1进行比较.

解 下面分别用分段二阶、三阶拉格朗日插值多项式和三次样条插值方法分别构造被积函数,计算数值积分.

方法1  (1) 用二阶拉格朗日插值多项式构造被积函数,输入程序

>>  X =[0,0.1571, 0.4712,0.5498,0.6283,0.8639, 1.0210,1.0996,1.5708];

Y=[0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000];

L1=lfun(X(1:3),Y(1:3)), L2=lfun(X(3:5),Y(3:5)),

 L3=lfun(X(5:7),Y(5:7)), L4=lfun(X(7:9),Y(7:9)),

运行后输出多项式L1,L2,L3,L4,然后输入程序

>> syms x

L1=inline('-14644281526214207/140737488355328000*x.^2+2280009553133670687/2251799813685248000*x');

L2=inline('-88810055001957943/351843720888320000*x.^2+15892292661979563/14073748835532800*x-15511422992738729/703687441776640000');

L3=inline('-261101353037957127/703687441776640000*x.^2+362054074498830351/281474976710656000*x-207985426949036897/2814749767106560000');

L4=inline('-131688789655946847/281474976710656000*x.^2+104193434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[Q141,FCNT141] = quad(L1,0, 0.4712,1.e-14);

[Q142,FCNT142] =quad(L2, 0.4712, 0.6283,1.e-14);

[Q143,FCNT143] = quad(L3, 0.6283, 1.0210,1.e-14);

[Q144,FCNT144] = quad(L4, 1.0210, pi/2,1.e-14);

Q14= Q141+ Q142+ Q143+ Q144,

FCNT14= FCNT141+ FCNT142+ FCNT143+ FCNT144

[Q41,FCNT41] = quad(L1,0, 0.4712,1.e-4);

[Q42,FCNT42] =quad(L2, 0.4712, 0.6283,1.e-4);

[Q43,FCNT43] = quad(L3, 0.6283, 1.0210,1.e-4);

[Q44,FCNT44] = quad(L4, 1.0210, pi/2,1.e-4);

Q4= Q41+ Q42+ Q43+ Q44,

FCNT4= FCNT41+ FCNT42+ FCNT43+ FCNT44

[Ql41,FCNTl41] = quadl(L1,0, 0.4712,1.e-4);

[Ql42,FCNTl42] =quadl(L2, 0.4712, 0.6283,1.e-4);

[Ql43,FCNTl43] = quadl(L3, 0.6283, 1.0210,1.e-4);

[Ql44,FCNTl44] = quadl(L4, 1.0210, pi/2,1.e-4);

Ql4= Ql41+ Ql42+ Ql43+ Ql44,

FCNTl4= FCNTl41+ FCNTl42+ FCNTl43+ FCNTl44

[Ql141,FCNTl141] = quadl(L1,0, 0.4712,1.e-14,2);

[Ql142,FCNTl142] =quadl(L2, 0.4712, 0.6283,1.e-14,2);

[Ql143,FCNTl143] = quadl(L3, 0.6283, 1.0210,1.e-14,2);

[Ql144,FCNTl144] = quadl(L4, 1.0210, pi/2,1.e-14,2);

Ql14= Ql141+ Ql142+ Ql143+ Ql144,

FCNTl14= FCNTl141+ FCNTl142+ FCNTl143+ FCNTl144

wQ14=double(abs(1- Q14)), wQ4=double(abs(1- Q4))

wQl4=double(abs(1- Ql4)), wQl14=double(abs(1- Ql14))

运行后屏幕显示如下

Q14 =                FCNT14 =     Q4 =                FCNT4 =

0.99957595925413        52       0.99957595925413          52

Ql4 =                FCNTl4 =      Ql14 =              FCNTl14 =

0.99957595925413        72       0.99957595925413          72

wQ14 =                             wQ4 =

4.240407458748763e-004             4.240407458748763e-004

wQl4 =                              wQl14 =

4.240407458748763e-004              4.240407458748763e-004

方法2  用三阶拉格朗日插值多项式构造被积函数,输入程序

>>  X =[0,0.1571, 0.4712,0.5498,0.6283,0.8639, 1.0210,1.0996,1.5708];

Y=[0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000];

L1=lfun(X(1:4),Y(1:4)), L2=lfun(X(4:7),Y(4:7)),

 L3=lfun(X(7:9),Y(7:9)),运行后输出多项式L1,L2,L3.然后输入程序:

>> syms x

L1=inline('-9070382563990323/56294995342131200*x.^3-793956227651389/281474976710656000*x.^2+2253151961697736691/2251799813685248000*x');

L2=inline('-41118170146065463/351843720888320000*x.^3-1700780706742359/21990232555520000*x.^2+735208215754549839/703687441776640000*x-25688208496779769/2814749767106560000');

L3=inline('-131688789655946847/281474976710656000*x.^2+104193434945799191/70368744177664000*x-9309740592764125768759/54295819320043765760000');

[Q141,FCNT141] = quad(L1,0, 0.5498,1.e-14);

[Q142,FCNT142] =quad(L2, 0.5498, 1.0210,1.e-14);

[Q143,FCNT143] = quad(L3, 1.0210, pi/2,1.e-14);

Q14= Q141+ Q142+ Q143 ,FCNT14= FCNT141+ FCNT142+ FCNT143

[Q41,FCNT41] = quad(L1,0, 0.5498,1.e-4);

[Q42,FCNT42] =quad(L2, 0.5498, 1.0210,1.e-4);

[Q43,FCNT43] = quad(L3, 1.0210, pi/2,1.e-4);

Q4= Q41+ Q42+ Q43, FCNT4= FCNT41+ FCNT42+ FCNT43

[Ql41,FCNTl41] = quadl(L1,0, 0.5498,1.e-14);

[Ql42,FCNTl42] =quadl(L2, 0.5498, 1.0210,1.e-14);

[Ql43,FCNTl43] = quadl(L3, 1.0210, pi/2,1.e-4);

Ql4= Ql41+ Ql42+ Ql43, FCNTl4= FCNTl41+ FCNTl42+ FCNTl43

[Ql141,FCNTl141] = quadl(L1,0, 0.5498,1.e-14);

[Ql142,FCNTl142] =quadl(L2, 0.5498, 1.0210,1.e-14);

[Ql143,FCNTl143] = quadl(L3, 1.0210, pi/2,1.e-4);

Ql14= Ql141+ Ql142+ Ql143, FCNTl14= FCNTl141+ FCNTl142+ FCNTl143

wQ14=double(abs(1- Q14)), wQ4=double(abs(1- Q4))

wQl4=double(abs(1- Ql4)), wQl14=double(abs(1- Ql14))

运行后屏幕显示如下

Q14 =                 FCNT14 =   Q4 =                 FCNT4 =

0.99975250938584       39       0.99975250938584          39

Ql4 =                  FCNTl4 =   Ql14 =               FCNTl14 =

0.99975250938584        54      0.99975250938584          54

wQ14 =                            wQ4 =

2.474906141644029e-004            2.474906141644029e-004

wQl4 =                             wQl14 =

2.474906141642919e-004            2.474906141642919e-004

方法3  三次样条插值多项式计算数值积分.

输入程序

>> a=0; b= pi/2;

X=[0,0.1571, 0.4712,0.5498,0.6283,0.8639, 1.0210,1.0996,1.5708];

fs=[0,0.1565,0.4540,0.5225,0.5878,0.7604,0.8526,0.8910,1.0000];

pp=spline(X,fs);

Qsl=quadl(@ppval,a,b,[],[],pp),wQs1=double(abs(1-Qsl))

Qs8=quad8(@ppval,a,b,[],[],pp),wQs8=double(abs(1-Qs8))

Qs=quad(@ppval,a,b,[],[],pp),wQs=double(abs(1-Qs))

运行后屏幕显示结果如下

Qsl =                      Qs8 =                 Qs =

1.00012614786269        1.00012598484705    1.00012589347894

wQs1 =                     wQs8 =                  wQs =

1.261478626870538e-004 1.259848470498426e-004  1.258934789407640e-004

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