矩陣與數(shù)值分析

本文格式為Word版，下載可任意編輯——矩陣與數(shù)值分析

矩陣與數(shù)值分析

學(xué)院專業(yè)班級(jí)學(xué)號(hào)姓名

電子信息與電氣工程學(xué)部

生物醫(yī)學(xué)工程

劉江濤

1：考慮計(jì)算給定向量的范數(shù)；輸入向量x?(x1,x2,?,xn)T，輸出x1,x2,x?，請(qǐng)編制一個(gè)通用程序，并用你編制的程序計(jì)算如下向量的范數(shù)：

1??11Tx??1,,,?,?,y??1,2,?,n?

n??23對(duì)n?10,100,1000甚至更大的n計(jì)算其范數(shù)，你會(huì)發(fā)現(xiàn)什么結(jié)果？你能否修改你的程序使得計(jì)算結(jié)果相對(duì)確切呢？

通用求范數(shù)程序：functionNORM(x)y1=sum(abs(x));y2=(sum(x.^2))^(1/2);y3=max(abs(x));

fprintf('1-范數(shù)=%g；2-范數(shù)=%g；inf-范數(shù)=%g\\n',y1,y2,y3);例題的運(yùn)行程序：functionxianglaing(n)x=[];y=[];fori=1:nx(i)=1/i;y(i)=i;end

disp('x的范數(shù)：');NORM(x');disp('')

disp('y的范數(shù)：');NORM(y');運(yùn)行結(jié)果如下表：

Tn范數(shù)101001000x12.928975.187387.48547x2x?1111.24491.278661.28216n范數(shù)101001000y1555050500500y2y?10100100019.6214581.67918271.1根據(jù)上述的兩個(gè)表的運(yùn)行結(jié)果，我們可以得知無(wú)論n的值如何變化，對(duì)于x??1恒成立；y??n恒成立，其1-范數(shù)與2-范數(shù)隨著n的增大而增大，但是其變化越來(lái)越小，這是由于計(jì)算在進(jìn)行數(shù)值計(jì)算時(shí)有誤差存在，對(duì)于表達(dá)式(1)當(dāng)n很大時(shí)

1卻很n小，會(huì)出現(xiàn)“大數(shù)吃小數(shù)的現(xiàn)象〞；修改方案：當(dāng)n很大時(shí)我們避免用n做除數(shù)，由于當(dāng)n十分大時(shí)

1?0成立；所以在求解其范數(shù)時(shí)我們從小數(shù)開始相加，無(wú)窮個(gè)十分n小的數(shù)值相加也可能是個(gè)很大的數(shù)，從而可以避免兩個(gè)數(shù)相加時(shí)出現(xiàn)“大數(shù)吃小數(shù)〞的現(xiàn)象；

2：考慮y?f(x)?ln(1?x)，其中定義f(0)?1，此時(shí)f(x)是連續(xù)函數(shù)，用此公x式計(jì)算當(dāng)x?[?10?15,10?15]時(shí)的函數(shù)值，畫出圖像。另一方面，考慮下面算法：

d?1?x;ifd?1theny?1elsey?lnd/(d?1)endif

用此算法計(jì)算x?[?10?15,10?15]時(shí)的函數(shù)值，畫出圖像，比較一下發(fā)生了什么？程序：

x=-10^(-15):10^(-20):10^(-15);

if(x==0)f=1;elsef=log(1+x)/x;endfigure(1)plot(x,f);d=1+x;ifd==1y=1;else

y=log(d)/(d-1);endfigure(2)Plot(x,y);

有圖可知，直接用公式f(x)?ln(1?x)計(jì)算x?[?10?15,10?15]的函數(shù)值時(shí)，除了在xx?0出的值為1，其他的值都是無(wú)限趨近于1；而利用算法二算出的結(jié)果全為1；出

現(xiàn)這這狀況的原因是x的取值十分接近于0，在用公式d?1?x求d得過(guò)程中出現(xiàn)了大數(shù)吃小數(shù)的狀況，所以在用計(jì)算機(jī)計(jì)算時(shí)d?1恒成立，從而使y?1恒成立；

3：首先編寫一個(gè)利用秦九韶算法計(jì)算一個(gè)多項(xiàng)式在定點(diǎn)的函數(shù)值的通用程序，你的程序包括輸入多項(xiàng)式的系數(shù)以及定點(diǎn)，輸出函數(shù)值，利用你編寫的程序計(jì)算

f(x)?(x?2)9?x9?18x8?144x7?672x6?2023x5?4032x4?5376x3?4608x2?2304x?512在x?2鄰域附近的值，畫出p(x)在x?[1.95,20.5]上地圖像。

秦九韶算法的通用程序：

%A為多項(xiàng)式的以升冪排列的系數(shù)，x為初始值functionp=qinjiushao(A,x)a=A;[~,n]=size(a);n=n-1;S=[];

S(n+1)=a(n+1);fork=n:-1:1

S(k)=x.*S(k+1)+a(k);

endp=S(1);

利用上述程序計(jì)算p(x)在x?2鄰域附近的值具體見下表：

xp(x)1.95-9.6634e-131.96-1.1255e-111.972.8422e-121.987.3896e-121.995.3433e-122.000xp(x)2.01-3.7517e-122.02-7.7875e-122.034.4338e-122.04-3.6380e-122.055.1159e-12當(dāng)x?[1.95,20.5]時(shí)，p(x)的圖像如下：畫圖程序如下：functionhuatu(A,x)[~,n]=size(x);fori=1:n

y(i)=qinjiushao(A,x(i));endplot(x,y);程序運(yùn)行如下：>>x=1.95:0.01:20.5;

>>A=[-5122304-46085376-40322023-672144-181];>>huatu(A,x)

3x10112.521.5p(x)10.50-0.50510x152025

4：編制計(jì)算機(jī)給定矩陣A的LU分解和PLU分解的通用程序，然后利用你編寫的程序完成下面兩個(gè)計(jì)算任務(wù)：

考慮

?1??1?A??????1???10?1?????01??Rn?n

?11??11??0??????1??1自己取定x?Rn，并計(jì)算b?Ax。然后用你編制的不選主元的Gauss消去法求解

?。對(duì)n從5到30估計(jì)計(jì)算解的精度。該方程組，記你計(jì)算出的解為x（2）對(duì)n從5到30計(jì)算出其逆矩陣。LU分解的通用程序：functionLU(A)[m,n]=size(A);L=zeros(m,n);U=zeros(m,n);forj=1:n

U(1,j)=A(1,j);L(j,j)=1;

endforj=2:m

L(j,1)=A(j,1)/U(1,1);endfori=2:nforj=i:nsum=0;fork=1:i-1

sum=sum+L(i,k)*U(k,j);end

U(i,j)=A(i,j)-sum;endforj=i+1:nsum=0;fork=1:i-1

sum=sum+L(j,k)*U(k,i);end

L(j,i)=(A(j,i)-sum)/U(i,i);endenddisp('L=');disp(L);disp('U=');disp(U);

PLU分解的通用程序：functionPLU(A)

[~,n]=size(A);Ip=1:n;fork=1:n-1

[~,r]=max(abs(A(k:n,k)));r=r+(k-1);ifr>k

A([k,r],:)=A([r,k],:);Ip([k,r])=Ip([r,k]);endforp=k+1:nmu=A(p,k)/A(k,k);A(p,k)=mu;

A(p,k+1:n)=A(p,k+1:n)-mu*A(k,k+1:n);endendp=eye(n,n);P=zeros(n,n);fori=1:n

P(i,1:n)=p(Ip(i),1:n);end

L=tril(A,-1)+eye(n);U=triu(A);disp('P=')disp(P);disp('L=')disp(L);

disp('U=');disp(U);

我選取b?[123?n]Tn?[5,30],n?N*，則我們可以計(jì)算出其確切解

2n?1?12n?2?1?1x?[?n?1?n?2?2222n?1T]2n?1n?[5,30]n?N*

實(shí)現(xiàn)求解的程序如下：functionINV(n)A=ones(n);fori=2:nforj=1:i-1

A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendfori=1:nb(i)=i;

X(i)=-(2^(n-i)-1)/(2^(n-i));end

X(n)=(2^n-1)/(2^(n-1));[L,U]=LU(A);x1=inv(U)*inv(L)*b'disp(x1);

我們利用上述的程序計(jì)算出的結(jié)果如下表：

nx1,x2,x3,?,xn

567…-0.9375,-0.8750,-0.7500,-0.5000,1.9375-0.9688,-0.9375,-0.8750,-0.7500,-0.5000,1.9688-0.9844,-0.9688,-0.9375,-0.8750,-0.7500,-0.5000,1.9844…利用我們求出的確切解與用程序求出的近似解求其誤差，再利用matlab編程實(shí)現(xiàn)時(shí)其誤差為零。（2）

對(duì)于題目中的A的求逆的通用程序：functionINV(n)A=ones(n);fori=2:nforj=1:i-1

A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendB=inv(A);

disp('A的逆為：');disp(B);n=5時(shí)：

-0.2500-0.1250-0.0625-0.0625??0.5000??00.5000-0.2500-0.1250-0.1250???1?A??000.5000-0.2500-0.2500??0000.5000-0.5000????0.25000.12500.06250.0625?0.5000?n=6時(shí)：

-0.2500-0.1250-0.0625-0.0313-0.0313??0.5000?0?0.5000-0.2500-0.1250-0.0625-0.0625???0?00.5000-0.2500-0.1250-0.1250A?1???

000.5000-0.2500-0.2500?0??0?0000.5000-0.5000??0.50000.25000.12500.06250.03130.0313????n=7時(shí)：

-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156??0.5000?0?0.5000-0.2500-0.1250-0.0625-0.0313-0.0313???0?00.5000-0.2500-0.1250-0.0625-0.0625??A?1??0000.5000-0.2500-0.1250-0.1250?

?00000.5000-0.2500-0.2500???000000.5000-0.5000???0.50000.25000.12500.06250.03130.01560.0156???n=8時(shí)：

-0.2500-0.1250-0.0625-0.0313-0.0156-0.0078-0.0078??0.5000?0?0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156???0?00.5000-0.2500-0.1250-0.0625-0.0313-0.0313??0000.5000-0.2500-0.1250-0.0625-0.0625?A?1???00000.5000-0.2500-0.1250-0.1250???000000.5000-0.2500-0.2500???0000000.5000-0.5000???0.25000.12500.06250.03130.01560.00780.0078???0.5000?

在此不再一一列舉，都可以用上述程序算出；

5：編制計(jì)算對(duì)稱正定陣的Cholesky分解的通用程序，并利用你編制的程序計(jì)算

Ax?b，其中A?(aij)?Rn?n,aij?1，b可以有你自己取定，對(duì)n從10到20驗(yàn)

i?j?1證程序的可靠性。

Cholesky分解求L的通用程序：%LT代表L的轉(zhuǎn)置function[L,LT]=Cholesky(A)[n,m]=size(A);L=zeros(n,n);forj=1:msum=0;fork=1:j-1

sum=sum+(L(j,k))^2;end

L(j,j)=(A(j,j)-sum)^(1/2);fori=j+1:nsum=0;fork=1:j-1

sum=sum+L(i,k)*L(j,k);end

L(i,j)=(A(i,j)-sum)/L(j,j);endendL=L;LT=L';

求解Ax?b的方程解的程序如下：functioncholesky_qiu_jie(n,b)A=zeros(n);fori=1:nforj=1:n

A(i,j)=1/(i+j-1);endend

[L,LT]=Cholesky(A);Y=inv(LT)*inv(L)*b;disp(Y');>>n=10;

>>b=[1000000000]';>>cholesky_qiu_jie(n,b)1.0e+06*

Columns1through5

0.000099995224850

-0.004949584178747

0.079191090719606

-0.6005186308358132.522130446082295

Columns6through10

-6.305225874989181

9.607833253959356

-8.749889022961289

4.374900125373453-0.923581794866949

>>n=11;

>>b=[10000000000]';>>cholesky_qiu_jie(n,b)1.0e+07*Columns1through5

0.000012085421839-0.0007244602419490.014116768083723

-0.1316795188512550.690983916160226

Columns6through10-2.210252357004733

4.471585077358848

-5.747468746188864

4.548898816586863-2.021271620300206

Column11

0.385801129112326n=12;

>>b=[111111111111]';>>cholesky_qiu_jie(n,b)1.0e+08*Columns1through5-0.000000106586186

0.000015496760612

-0.000549242897797

0.008320236358643-0.067089571289062

Columns6through100.321428166562500

-0.969534469375000

1.888382749062500

-2.3698240640625001.849525708125000

Columns11through12

-0.8162375431250000.155564209082031n=13;

>>b=[1111111111111]';>>cholesky_qiu_jie(n,b)1.0e+10*Columns1through50.000000017255662

-0.000002801898660

0.000111021739307

-0.0018912192593750.017318444275000

Columns6through10

-0.0955644689000000.338540141400000-0.795937993200000

1.255304630000000-1.312819645600000

Columns11through13

0.873192546400000-0.3343542442000000.056103628050000>>n=14;

>>b=[11111111111110]';>>cholesky_qiu_jie(n,b)1.0e+15*Columns1through5-0.000000018973466

0.000003062445902

-0.000120619605275

0.002041840742992-0.018571410648157

Columns6through100.101712531234567

-0.357252637738148

0.831474891793232

-1.2948808998775841.331275372322000

Columns11through14-0.862636505012200-0.002171662407011

>>n=15;

>>b=[111111111111101]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through50.000000009159679

-0.000000881206821

0.000013799271738

0.314323195943382

-0.045197162809621

0.000103167625626-0.003987948970647

Columns6through100.038905667203020

-0.198056922710865

0.601726021228474

-1.1199202319749541.191735365493764

Columns11through15-0.450794442163960

-0.496191645217859

0.756621451313167

-0.3998383557115620.079684802839153

>>n=16;

>>b=[1111111111111010]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through50.000000106506398

-0.000014879116455

0.000507146024044

-0.0072578714099440.052560866946652

Columns6through10-0.198284962679135

0.290609837891524

0.581225032778005

-3.4971342401756436.758520234521732

Columns11through15-5.326111865614911

-2.146668037387352

8.700448141028517

-8.1619190411434813.588648336729691

Column16

-0.635128701903556>>n=17;

>>b=[11111111111110101]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through50.000000014777084

-0.000002118394303

0.000074232548837

-0.0010976069501010.008314003699029

Columns6through10-0.034029369522343

0.065917257262332

0.015959536689882

-0.3776216029047730.819475115416494

Columns11through15-0.660163184213223

-0.295786673257216

1.100782605917994

-0.9529071409935970.339855611203294

Columns16through17

-0.014551513983896-0.014219174116106>>n=18;

>>b=[111111111111101012]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through50.000000155967115

-0.000021645380087

0.000727005315139

-0.0101303914094080.069910333639844

Columns6through10-0.237692894985775

0.207119320267744

1.320710222270083

-5.3628332770859458.761578864695567

Columns11through15-5.416385894431453

-3.678202876264916

9.162910938126201

-7.6100281754920474.801505077240169

Columns16through18

-3.4602807804032691.887745959501761-0.436632060441746>>n=19;

>>b=[1111111111111010120]';>>cholesky_qiu_jie(n,b)1.0e+17*

Columns1through50.000000000186975

0.000000487130031

-0.000037337486618

0.000907963809187-0.010389239660544

Columns6through100.065117610127561

-0.236770662992548

0.490423724491871

-0.485761304739977-0.006945225699907

Columns11through150.375039361697087

-0.191724740660743

0.584539934693850

-1.6646332685136670.883570199543987

Columns16through191.944997967101147-0.384572610629205

n=20;

b=[11111111111110101200]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5-0.000000005297720

0.000002202326696

-0.000135104609897

-3.195222688027543

1.831459798885819

0.003035236689788-0.033417756031216

Columns6through100.204643883440640

-0.727554702023794

1.433949931025033

-1.101616652841183-1.111630914232476

Columns11through152.722140887157841

-0.945117835365089

-0.676679844961023

-1.6436973612162932.211735585037371

Columns16through203.295476378249488

-7.583406568516137

5.452181213933089

-1.6592423652438400.159333719326849

A?1利用cholesky分解求出的解的確切度高于直接，由于當(dāng)n逐漸增大的過(guò)程中，

越來(lái)越接近奇異矩陣，使得計(jì)算結(jié)果的誤差增大，而使用cholesky分解可以避免這種現(xiàn)象的產(chǎn)生，是計(jì)算結(jié)果更加確切。

6：（1）編制程序House(x)，其作用是對(duì)輸入的向量x，輸出單位向量u使得

(I?2uuT)x?x2e1。

編制Householder變換陣H?I?2uuT?Rn?n乘以A?Rn?m的程序HA，注意，你的程序并不顯式的計(jì)算出H。

考慮矩陣

?1???1A???2???10?0?232224??23?e??

??37?75/2??3用你編制的程序計(jì)算H使得HA的第一列為?e1，并將HA的結(jié)果顯示出來(lái)。編制House(x)，其作用是對(duì)輸入的向量x,輸出單位向量U使得

(I?2uuT)x?x2e1：

House(x)通用程序：functionHouse(x)[n,~]=size(x);e1=zeros(n,1);e1(1)=1;w=x-norm(x)*e1;U=w/sqrt((w'*w));disp('U=');disp(U);

編制Householder變換陣H?I?2uuT?Rn?n乘以A?Rn?n的程序HA，注意，你的

程序并不顯示的計(jì)算出H。

Householder變換陣H通用程序：functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);考慮矩陣

?1??-1A??-2???10?0?24??323?2e??

?2?37?2752??3用你編制的程序計(jì)算H使得HA的第一列為?e1的形式，并將HA的結(jié)果顯示：

HA的結(jié)果顯示的通用程序：

functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);

w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);程序運(yùn)行結(jié)果：

>>A=[1234;-13sqrt(2)sqrt(3);-22exp(1)pi;-sqrt(10)2-37;0275/2];>>Householder(A)HA=

4.0000-2.83111.4090-6.53780.00001.38960.8839-1.78050.0000-1.22081.6576-3.88360.0000-3.0925-4.6770-4.107802.00007.00002.5000

7：用jacobi和Gauss-Scidel迭代求解下面的方程組，輸出每一步的誤差xk?x*；

?5x1?x2?3x3??2???x1?2x2?4x3?1??3x?4x?15x?10123?Jacobi迭代通用程序：functionjacobi(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin=epfprintf('%d步誤差\\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;n=n+1;end

disp('近似解')disp(y)

在命令窗口輸入以下命令就可以輸出每一步的誤差及最優(yōu)近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;>>jacobi(A,b,x0,ep)

近似解

-0.081967231207363-1.8032786473803261.131147556193974Gauss-Scidel迭代法通用程序：functionGaussseidel(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin=epfprintf('%d步誤差\\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;

n=n+1;end

disp('x的近似最優(yōu)解：')disp(y);

在命令窗口輸入以下命令就可以輸出每一步的誤差及最優(yōu)近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;

>>Gaussseidel(A,b,x0,ep)x的近似最優(yōu)解：-0.081967202756916-1.8032785851349111.131147515484593

8：取不同的初值用Newton迭代法以及弦截法求解x3?2x2?10x?100?0的實(shí)根，列表或者畫圖說(shuō)明收斂性；

Newton迭代法通用程序：

%使用說(shuō)明f為符號(hào)表達(dá)式，x0為初始值functionNewton(f,x0)symsx;x=x0;

x1=x0-eval(f/diff(f));while(norm(x1-x0)>0.000001)x=x1;

xk=x1-eval(f/diff(f));x0=x1;

x1=xk;disp(xk);enddisp('xk=')disp(xk)

弦截法通用程序：

%使用說(shuō)明f為符號(hào)表達(dá)式，x0,x1為初始值functionXuanjiefa(f,x0,x1)symsx;x=x0;f1=eval(f);x=x1;f2=eval(f);

x2=x1-(f2/(f2-f1))*(x1-x0);while(norm(x2-x1)>0.000001)x=x1;f1=eval(f);x=x2;f2=eval(f);

xk=x2-(f2/(f2-f1))*(x2-x1);x1=x2;x2=xk;enddisp('xk=')disp(xk)

Newton法：

k0123xk33.51023.46113.4606f(xk)-252.99620.03013.1407e-06xk?1?xk0.51020.04915.0362e-04k01234xk24.13333.54043.46193.4606f(xk)-6446.11784.85320.077442.0766e-5xk?1?xk2.13330.59290.07860.0013k0123xk43.51363.46123.4606f(xk)363.19820.03424.059e-6xk?1?xk0.48650.05245.726e-4弦截法：

k012345xk33.23.48793.45913.46063.406f(xk)-25-14.75201.6418-0.0907-5.1186e-041.6105e-07xk?1?xk0.20.28790.02880.00158.5667e-6k0123456xk233.64103.44283.45993.46063.4606f(xk)-64-2511.1937-1.0607-0.03901.4420e-4-1.9486e-8xk?1?xk10.64100.19830.01726.5469e-62.4129e-6從上述兩個(gè)表的結(jié)果我們可以對(duì)照得出，在不同的初始值的條件下運(yùn)行程序都具有收斂性，但是不同的初始值有不同的收斂速度；

9：用二分法求解方程excosx?2?0在區(qū)間[0,4?]上所有根。functioner_fen_fa(f,a,b)symsx;

e=0.0000001;while((b-a)>e)m=(b+a)/2;x=a;f1=eval(f);x=m;f2=eval(f);iff1*f2>=0a=m;elseb=m;endendc=(a+m)/2;disp(c);

根據(jù)函數(shù)的特點(diǎn)將其定義域分為4等分分別為[0,?]，[?,2?]，[2?,3?]，[3?,4?]分別運(yùn)行程序：>>symsx

>>f=exp(x)*cos(x)+2;>>a=0;>>b=pi;>>er_fen_fa(f,a,b近似零點(diǎn)：1.8807>>symsx

>>f=exp(x)*cos(x)+2;

a=pi;>>b=2*pi;>>er_fen_fa(f,a,b)近似零點(diǎn)：4.6941>>symsx

>>f=exp(x)*cos(x)+2;>>a=2*pi;>>b=3*pi;>>er_fen_fa(f,a,b)近似零點(diǎn)：7.8548>>symsx

>>f=exp(x)*cos(x)+2;>>a=3*pi;>>b=4*pi;>>er_fen_fa(f,a,b)近似零點(diǎn)：10.9955

;x2?4.6941;x2?7.8548;x4?10.9955由上述程序可知其全部解為x1?1.880710：考慮函數(shù)f(x)?sin(?x),x?[0,1]。用等距節(jié)點(diǎn)作f(x)的Newton插值，求出插值多項(xiàng)式以及f(x)的圖像，觀測(cè)收斂性。

Newton插值求插值多項(xiàng)式的通用程序：functionP=Nweton_xhazhi(X,f)symsx;p=0;

[~,n]=size(X);fori=1:n-1F=0;w=1;sum=1;forl=1:i+1;w=w*(x-X(l));endforj=1:i+1x=X(j);

F=F+f(j)/(eval(diff(w)));endfort=1:isymsx;

sum=sum*(x-X(t));endp=p+F*sum;endP=p+f(1);Disp(‘P=’);

用等距節(jié)點(diǎn)作f(x)的Newton插值多項(xiàng)式的通用程序：X=0:0.1:1;f=sin(pi*X);

O=Nweton_xhazhi(X,f);Y=simple(O);disp(Y);

plot(X,f)P=

x*(11153755408641612500000*x^919582893765599948750000*x^76457490747185549971250*x^5452380176499959382000*x^3

++-

-

55768777028538456250000*x^8256281087074564729875000*x^61150585140486343857745625*x^42327269296084898261641325*x^2

++++

4715634276489610818*x-1414847701376318005254584)

畫圖程序：X=0:0.1:1;x1=0:0.0001:1;f1=sin(x1*pi);symsxx0=0:0.001:1;f=sin(pi*X);

O=Nweton_xhazhi(X,f);x=x0;Y=eval(O);plot(x0,Y,'r');holdonplot(x1,f1);

1.210.80.60.40.20-0.200.10.20.30.40.50.60.70.80.91

由上圖可知：紅色圖形是利用插值多項(xiàng)式畫出的，藍(lán)色圖像是利用原函數(shù)畫出的，將兩者對(duì)比我們可以得出結(jié)論：該插值多項(xiàng)式具有收斂性

11：對(duì)函數(shù)f(x)?1,x?[?5,5]，取不同的節(jié)點(diǎn)數(shù)n,用等距節(jié)點(diǎn)作Lagrange1?x2插值，觀測(cè)Runge現(xiàn)象。

Lagrange插值通用程序：

%X為節(jié)點(diǎn)，y為節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值，x0為插值節(jié)點(diǎn)functionLagrange(X,y,x0)symsx;[~,n]=size(X);sum=0;fori=1:n

sum=sum+y(i)*L(X,i);endx=x0;p=eval(sum);plot(x0,p,'r-')holdon

plot(X,y,'b-');

functionl=L(X,k)symsx;sum1=1;sum2=1;[~,n]=size(X);if(k==1)fori=k+1:n

sum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));l=sum1/sum2;endelse

fori=1:k-1

sum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));endfori=k+1:n

sum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));end

l=sum1/sum2;end

取不同的節(jié)點(diǎn)數(shù)n,用等距節(jié)點(diǎn)作Lagrange插值，觀測(cè)Runge現(xiàn)象:

當(dāng)n=10;

X=-5:0.1:5;f=1./(1+X.^2);x0=-5:5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

當(dāng)n=20;

X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

當(dāng)n=50;

X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.2:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

n=100;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.1:5;

Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

由以上幾個(gè)圖形我們可以得知，插值節(jié)點(diǎn)越多時(shí)，插值多項(xiàng)式的收斂性越好；12：令f(x)?e3xcos(?x)，考慮積分?2?0f(x)dx。區(qū)間分為50，100，200，500，

1000等，分別用復(fù)合梯形以及復(fù)合Simpson積分公式計(jì)算積分值，將數(shù)值積分的結(jié)果與確切值比較，列表說(shuō)明誤差的收斂性。

復(fù)合梯形公式通用程序：

%f為被積函數(shù)，n為積分區(qū)間等分的數(shù)目，a,b分別為積分下上限functiontixing(f,n,a,b)symsx;sum=0;xk=a:(b-a)/n:b;fori=2:n-1x=xk(i);

sum=sum+eval(f);symsx;endx=a;f1=eval(f);x=b;f2=eval(f);

T=(b-a)/(2*n)*(f1+f2+2*sum);disp(T)

復(fù)合sipson公式：functionSimpson(f,n,a,b)symsx;sum=0;sum1=0;xk=a:(b-a)/n:b;fori=1:n-1

x=(xk(i)+xk(i+1))/2;sum=sum+eval(f);symsx;endfori=2:n-1x=xk(i);

sum1=sum1+eval(f);symsx;endx=a;

f1=eval(f);x=b;f2=eval(f);

T=(b-a)/(6*n)*(f1+f2+2*sum1+4*sum);disp(T)

利用matlab內(nèi)置函數(shù)int()求出該積分式的真是值為：>>symsx

>>int(exp(3*x)*cos(pi*x),0,2*pi)ans=

3.5232e+07

利用符合simpson公式計(jì)算如下表：

n501002005001000200030005000近似值FK2.3147e+072.9066e+073.2156e+073.4010e+073.4623e+073.4928e+073.5030e+073.5111e+07誤差F-Fk1.2085e+076.1664e+063.0761e+061.2228e+066.0960e+053.0430e+052.0275e+051.2159e+05

利用復(fù)合梯形公式計(jì)算如下表：

n501002005001000200030005000近似值FK2.3478e+072.9054e+073.2139e+073.4005e+0