大连理工大学矩阵与数值分析上机作业13382

1、矩阵与数值分析上机作业大连理工大学学校:学院:班级:姓名,学号：授课老师:注：编程语言Matlab1.考虑计算给定向量的范数：输入向量上=(£】,七2.,上产,翰出肛|1,同2,同8 请嫡制一个通用程序,并用你编制的程序计算如下向量的范数：Xtn = 10, 100, 1000甚至更大的计算其范数,你会更现什么结果？你能否修改你的程序 使得计算结先相好耕确呢？Norm, m函数 function s=Norm(x,m) +求向量X的范数取1,2, inf分别 表示1,2,无穷范数 n=length(x);s=0;switch mcase 1 +1-范数for i=l:ns=s+abs

2、(x(i);endcase 2 +2-范数for i=l:ns=s+x(i)A2;ends=sqrt(s);case inf学无穷-范数s=max(abs(x);end计算向x, y的范数Testi.mclear all;clc;nl=10;n2=100;n3=1000;xl=l./l:nlf;x2=l./l:n21;x3=l./1:n31;yl=l:nl,;y2=l:n2,;y3=l:n31;disp (1 n=10 时,);dispx 的 1-范数：,);disp (Norm(xl, 1);disp (' x 的 2-范数：1 ) ; disp (Norm (xlz 2);disp

3、 ( 1 x 的无穷-范数：1) ;disp (Norm (xlr inf);disp( *y 的 1-范数：1) ;disp (Norm(ylz 1);disp(y 的 2-范数：1) ;disp(Norm(yl,2);disp (1 y 的无穷-范数：,) ;disp (Norm (ylr inf);disp(1n=100 时,);disp (1 x 的 1-范数：1) ;disp (Norm(x2,1);disp ("的 2-范数:1 ) ; disp (Norm (x2z 2);disp (1 x 的无穷一范数:1) ;disp (Norm(x2r inf);disp( 1y

4、 的 1 -范数：1 ) ;disp (Norm(y2z 1);disp(' y 的 2-范数:1 ) ;disp (Norm (y2f 2);disp (1 y 的无穷-范数：1) ;disp (Norm (y2r inf);disp(1n=1000 时')；disp (的 1-范数:1) ; disp (Norm (x3,1);disp (1 x 的 2-范数：1) ; disp (Norm (x3,2);dispx 的无穷-范数：1) ;disp (Norm (x3, inf);disp ( 1 y 的 1-范数：');disp (Norm(y3z 1);disp

5、( fy 的 2-范数：1 ) ;disp (Norm(y3,2);disp ( 1 y 的无穷-范数：,);disp (Norm (y3r inf);运行结果：n=10 时x的范数:2.9290； x的2-范数：1.2449； x的无穷-范数：1V的范数：55； v的2-范数：19. 6214； y的无穷-范数：10n=100 时x的1-范数5 1874； x的2-范数：1.2787； x的无穷-范数：1y的范数：5050；y的2-范数:581.6786； y的无穷-范数：10021000时x的范数7. 4855； x的2-范数：1.2822；x的无穷-范数:1y 的 1-范数：500500

6、； y 的 2-范数：1.8271e+004； y 的无穷-范数：10002.考虑y = /(M)= 吗且,其中定义0)= 1.此时/(工)是连续函数.用此公式计算 当工£-10-叫10-叫时的函数值,画出图像.另一方面,考虑下面算法：d = 1 + 1tfd= theny=lelsey = nd/(d 1)end if用此算法计算上10-15J时的函数值,画出困像.比拟一下发生了什么？Test2. mclear all;clc;n=100;%|xfHh=2«10八(-15) /n;+步长x=-L0A(-15):h:10A(-15);学第一种原函数fl=zeros(1,n+

7、1);for k=l:n+lif x(k)=0fl(k)=log(1+x(k)/x(k);elsefl(k)=l;卑微如螃蚁、坚强似大象endendsubplot(2,1,1);ploz(k,flf1-r1);axis(-10A(-15)r10A(-15),-1,2);legend (1 原图 1);+第二种算法f2=zeros(lr n+1);for k=l:n+1d=l+x(k);if(d=l)f2(k)=log(d)/(d-1);elsef2(k)=l;endendsubplot(2rlr 2);plot(x,f2,'-r');axis(-10A(-15),10A(-15

8、),-1,2);legend (,第二种克法,)；运行结果,显然第二种算法结果不准确,是由于计算机中的舍入误差造成的,当xw-10F,10,时,d = l+K, 计算机进行舍入造成d恒等于1,结果函数值恒为1.3.首先编制一个利用暴九铝算法计算一个多项式在给定点的函数值的通用程序,你 的程序包括输入多项式幽系数以及给定点,榆出函数值.利用你编制的程序计算P(T)= (z-2)9=i9- 18j:8 4 144jt7 - 672/ + 2021/ _1032/ + 5376/ _ 46O8j:2 + 23O4x-512在/= 2邻城附近的值.画出(工)在j： W 1,95,20,5上的图像.W*

9、：秦九韶算法:QinJS. mfunction y=QinJS(a,x)* y输出函数值* a多项式系数,由高次到零次* x给定点n=length(a);s=a(l);for i=2:ns=s*x+a(i);endy=s;计j|p(x) :test3. mclear all;clc;x=1.6:0.2:2.4;%x=2 的邻域disp(,x=2 的邻城：1) ;xa=l -18 144 -672 2021 -4032 5376 -4608 2304 -512;p=zeros(1/5);for i=l:5p(i)=QinJS(ar x(i);enddisp相应多项式p值:*) ;pxk=L.95

10、:0.01:20.5;nk=Length(xk);pk=zeros(1, nk);k=l;for k=l:nkpk(k)=QinJS(a, xk(k);endploz(xk,pk,1-r1);xlabel (fxf) ; ylabel (1 p (x) f);运行结果：x=2的邻域： x =1.60001.80002. 00002. 20002. 4000相应多项式P值：P =1.0e-003 *-0.2621-0.000500.00050.2621p(x)在 X WL95,20.5上的图像4.编制计算给定矩阵百的AU分解和PLU分解的通用程序,然后用你编制的程序完 成下面两个计算任务：(1)

11、考虑Toof- 1 : 4= 一.Q1 £R"n- 1 -1 1 1- 1-1-11 自己取定上GR7并计算5= Ai.然后用你编制的不选主元和列主元的Caws消去法求解 该方程组,记你计算出的解为人对,从5到3.估计计算解的精度.(2)对n从5到和计算其逆矩阵.皿LU 分解,LUDecom. m function L,U=LUDecom(A) +不带列主元的LU分解 N = size(A);n = N (1);L=eye(n);U=zeros (n); for i=l:nu(l,i)=A(l,i) ;L(ir l)=A(i,l) /u(lzl); end for i=2:

12、nfor j=i:n z=0; for k=l:i-lz=z+L(irk)*U(kf j); endU(irj)=A(irj)-z;endfor j=i+l:nz=0;for k=l:i-lz=z+L (j , k) *U (k, i);endL(j,i) = (A(j,i)-z)/U(i,i);endendPLU 分解,PLUDecom. mfunction P,L,U =PLUDecom(A)%带列主元的LU分解mzm=size(A);U=A;P=eye(m);L=eye(m);for i=l:mfor j=i:mt(j)=U(j,i)；for k=l:i-lt (j) =t (j) -U

13、 (j , k) *U(k, i);endenda=i;b=abs(t(i);for j=i+l:mif b<abs(t(j)b=abs(t(j );a=j；endendif a=ifor j=l:mc=U(i,j);U(i,j)=U(a,j);U(azj)=c;endfor j=l:mc=P(ifj);P(i,j)=P(a,j);P(azj)=c;endc=t(a);t(a)=t (i);t(i)=c;endU (i, i) =t (i);for j=i+l:mu (j, i)=t (j ) /t (i);endfor j=i+l:mfor k=l:i-lU(i, j)=U(i,j)-

14、U(i,k)*U(k,j)；endendendL=tril(U,-1)+eye(m);U=triu(U,0);(1) (2)程序：Test4. mclear all;clc;for n=S:30x=zeros(n,1);A=-ones(n);A(:,n)=ones(n,1);for i=l:nA(i,i)=l；for j = (i+l) : (n-1)A(i,j)=0;endx(i)=l/i;enddisp ( 1 n=1 ) ; disp (n);disp方程精确解:;Xb=A*x; *系数 bdisp (,利用LU分解方程组的解:')；L, U =LUDecom (A) ; %LU

15、 分解xLU=U(Lb)disp (,利用PLU分解方程组的解:')；P,LfU =PLUDecom(A) ;%PLU 分解xPLU=U(L(Pb)+求解A的逆矩阵disp (,A的准确逆矩阵：,)；InvA=inv(A)InvAL=zeros (n) ; e利用LU分解求A的逆矩阵I=eye(n);for i=l:nInvAL(:ri)=U(LI(:fi);enddisp (,利用LU分解的A的逆矩阵：,)；InvALEnd运行结果：(1)只列出n=5, 6, 7的结果当n= 5方程精确解：卑微如蟋蚁、坚强似大象X =1.00000. 50000. 33330. 25000. 200

16、0利用LU分解方程组的解：xLU =1.00000. 50000. 33330. 25000. 2000利用PLU分解方程组的解:xPLU =1.00000. 50000. 33330. 25000. 2000当n=6方程精确解：x =1.00000. 50000. 33330. 25000. 20000. 1667利用LU分解方程组的解：xLU =1.00000. 50000. 33330. 25000. 20000. 1667利用PLU分解方程组的解:xPLU =1.00000. 50000. 33330. 25000. 20000. 1667当n= 7方程精确解：x =1.00000.

17、50000. 33330. 25000. 20000. 16670. 1429利用LU分解方程组的解：xLU =1.00000. 50000. 33330. 25000. 20000. 16670. 1429利用PLU分解方程组的解：xPLU =1.00000. 50000. 33330. 25000. 20000. 16670. 14292只列出n=5,6,7时A的逆矩阵的结果卑.微如蟋蚁、坚强似大象当n= 5A的准确逆矩阵:0. 5000-0. 2500-0. 1250-0.0625-0.06250. 5000-0. 2500-0. 1250-0.12500. 5000-0. 2500-0

18、.25000. 5000-0.50000. 50000. 25000. 12500. 06250.0625利用LU分解的A的逆矩阵:InvAL =0. 5000-0. 2500-0. 1250-0. 0625-0.062500. 5000-0. 2500-0. 1250-0.1250000. 5000-0. 2500-0.25000000. 5000-0.50000. 50000. 25000.12500. 06250.0625当n= 6A的准确逆矩阵:0. 5000-0. 2500-0. 1250-0.0625-0.03130. 5000-0. 2500-0. 1250-0.06250. 5

19、000-0. 2500-0.12500. 5000-0.25000.50000. 50000. 25000. 12500. 06250.0313-0. 0313-0.0625-0. 1250-0.2500-0.50000. 0313利用LU分解的A的逆矩阵:InvAL =0. 500000000. 5000-0. 25000. 50000000. 2500-0. 1250-0. 25000.5000000. 1250-0. 0625-0. 1250-0. 25000. 500000. 0625-0.0313-0.0625-0.1250-0.25000.50000.0313-0. 0313-0.

20、 0625-0. 1250-0. 2500-0. 50000. 0313当n= 7A的准确逆矩阵:InvA =0. 5000-0. 2500-0. 1250-0. 0625-0.0313-0.0156-0.015600. 5000-0.2500-0. 1250-0.0625-0. 0313-0. 0313000.5000-0. 2500-0.1250-0. 0625-0.06250000. 5000-0.2500-0.1250-0.125000000.5000-0.2500-0.2500000000. 5000-0.50000. 50000. 25000. 12500. 06250.03130

21、.01560.0156利用LU分解的A的逆矩阵:InvAL =-0. 0625-0.0313-0.0156-0.0156-0. 1250-0.0625-0. 0313-0. 0313-0. 2500-0.1250-0.0625-0.06250. 5000-0.2500-0. 1250-0.125000.5000-0. 2500-0.25000. 5000-0. 2500-0.125000.5000-0.2500000.5000000000000000.5000 -0.50000. 50000. 25000. 12500. 06250. 03130. 01560. 01565.编制计算时称正定阵

22、的6.及果y分解的通用程序,并用你绢制的程序计算Ar =4 其中,4 =(旬) Rnxn,旬=哥z?.b可以由你自己取走,对R从10到2U聆证程序的可靠性.皿Cholesky 分解：Choi esky. mfunction L=Cholesky(A)N = size(A);n = N (1);L=zeros (n);L(l,l)=sqrt(A(lrl);for i=2:nL(i,l)=A(i,l)/L(l,l)；endfor j=2:nsl=0;for k=l:j-1sl=sl+L(j,k)A2;endj)=sqrt(A(j, j)-sl);for i=j+l:ns2=0;for k=l:j-

23、1s2=s2+L(i,k)*L(j,k);endL(ir j) = (A(i, j)-s2) /L(j, j);endend计算 Ax=b； Test5. mclear all;clc;for n=10:20A=zeros(nr n);b=zeros (n,1);for i=l:nfor j=l:nA(izj)=l/(i+j-l);endb(i,l)=i；enddisp ( f n= 1 ) ; disp (n);disp (1方程组原始解,);xO=Abdisp (»利川Cholesky分解的方程组的解,);L=Cholesky(A)x=L、Lb)end运行结果：里微如螃蚁、坚强似

24、大象只列出了 n=10,11的结果n=10方程组原始解xO =1.0e+008 *-0.00000. 0010-0. 02330. 2330-1.21083. 5947-6. 32336.5114-3.62330. 8407利用Choi esky分解的方程组的解 x =1.0e+008 *- 0. 00000. 0010- 0. 02330. 2330-1.2105- 6. 32196.5100- 3. 62250. 8405n=11方程组原始解x0 =1.0e+009 *0. 0000-0.00020. 0046-0. 05670. 3687-1.40393. 2863-4.78694. 2

25、260-2. 06850. 4305利用Choi esky分解的方程组的解X =1.0e+009 *0. 0000-0. 00020. 0046-0.05630. 3668-1.39723. 2716-4. 76694. 2094-2. 06080. 42906. (1)编制程序月.!?0 其作用是对输入的向量上,£丁出单位向量“使得(/-2履二 同型.(2)编制班山店idder变换阵"=/一21/那'来以4股"小的程序4 注意,你的程 序并不显式的计算出.(3)考虑矩阵/1234 -135/2>/3A =22er ,一加2-37I0275/2)用你

26、编制的程序计算H使得HA的第一列为cei的形式,并将HA的结果显示.(1) House, mfunction u=House(x)n=length(x);el=eye(n,1);w=x-norm(xz 2)*el;u=w/norm(w,2);(2) Hou_A. mfunction HA=Hou_A(A)卑微如螃蚁、坚强似大象al=A(:f1);n=length (al);el=eye (n,1);w=al-norm(al/2)*el;u=w/norm(w,2);H=eye(n)-2*u*u1HA=H*A；test6.mclear all;clc;A=L 234;- 1 2 sqrt(2) s

27、qrt(3);- 2 2 exp (1) pi;- sqrt(10) 2 -3 7;027 5/2;HA=Hou_A (A)运行结果：H =0. 2500-0. 2500-0. 5000-0. 79060-0.25000.9167-0. 1667-0. 26350-0. 5000-0. 16670. 6667-0. 52700-0.7906-0.2635-0. 52700. 1667000001.00004. 0000-2. 58111.4090-6. 53780. 00000. 47300. 8839-1.78050. 0000-1.05411.6576-3. 88360. 0000-2.

28、8289-4. 6770-4. 10782. 00007. 00002. 50007.用Jaco儿和代求解下面的方程纨.输出迭代每一的误甚|上小 -上/|：(5X112 3X3 = 2一0+ 2上2 + 413 =1呢+ 42 + 15T3 = 10皿Jacobi 迭代：Jaccobi.function xr n=Jaccobi(A,b,xO)学-方程组系数阵A%方程组右端顶b%初始值xOv一求解要求精确度eps*一迭代步数限制M返回求得的解x*-返回迭代步数nM=1000;eps=l.Oe-5;D=diag (diag (A) ) ;*求A的对角矩阵L=-zril (A, -1) ; *求A

29、的下三角阵U=-triu (A, 1) ;*求A的上三角阵J=D(L+U);f=Db;卑微如螃蚁、坚强似大象x=J*xO+fn=l; +迭代次数err=norm(x-xO,inf)while(err>=eps)xO=x;x=J*xO+fn=n+l;err=norm(x-xOf inf)if(n>=M)disp (* Warning:迭代:多,可能不收敛T);return;endendGauss.SeideI 迭代：Gauss.Seidel.mfunction xf n=Gauss_Seidel(A,b,xO)%Gauss-Seidel迭代法解线性方程组+一方程组系数阵A*一方程组右

30、端项b%-初始值xO%一求解要求的精确度eps%迭代步数限制M返回求得的解X %-返回迭代步数n eps=l.Oe-5;M=10000;D=diag (diag (A) ) ; *求A的对角矩阵L=-tril (A, -1) ; *求A的下三角阵U=-riu (A, 1) ; %求A的上三角阵G= (D-L) U；f=(D-L)b；x=G*xO+fn=l;%迭代次数err=norm(x-xO,inf)while(err>=eps)xO=x;x=G*xO+fn=n+l;err=norm(x-xO,inf)if(n>=M)disp (1 Warning:迭代次数太多,可能不收敛！,)；

31、 return;endend解方程组,test7. mclear all;clc;A=(5 -1 -3;-1 2 4;-3 4 15;b=-2;l;10;disp ,精确解,;x=Abdisp迭代初始值,；xO=O;0;0disp ' Jacobi 迭代过程：1;xj,nj=JaccobiA,b,xO;disp ' Jacobi最终迭代结果：,； xjdisp ,迭代次数,；njdisp 1 Gauss-Seidel 迭代过程:1;xg,ng=Gauss_SeidelAzb,x0;disp 1 Gauss-Seidel 最汴迭代二 1 果：, ; xgdisp ,迭代次数,；

32、ng运行结果：精确解X =-0. 0820-1.80331.1311迭代初始值xO =000Jacobi迭代过程：x =-0. 40000. 50000. 6667err =0. 6667x =0. 1000-1.03330. 4533 err =1.5333 x =-0. 0820-1.80331.1311 err =9.6603e-006Jacobi 终迭代结果： xj =-0. 0820-1.80331.1311迭代次数nj =281Gauss-Sei do I 迭代过程: x =-0. 4000o. 3000p. 5067 err =q. 5067 x =-0.0360-0. 5313

33、0. 8012 err =0. 8313X =-0. 0256-1.11510. 9589 err =0. 5838 X =-0.0820-1.80331.1311 err =9.4021e-006Gauss-Seidel最终迭代结果： xg =-0. 0820-1.80331.1311 迭代次数 ng =208.取不同的初值用Nowkni迭代以及弦栈法求方程： + 2 + 10上_100=0的实根,列 表或者画图说明收敛速度.Newton 迭代法：Newton iter. mfunction xf iter,fvalue=Newtoniter(£,dfr xO,epsrmaxite

34、r) 卑微如螃蚁、坚强似大象*牛顿法X得到的近似解e iter迭代次数%fvalue函数在x处的值%f. df被求的非线性方程及导函数%xO初始值%eps允许误差限%maxiter最大迭代次数fvalue=subs(£,xO);dfvalue=subs(dff xO);for iter=l:maxiterx=xO-fvalue/dfvalueerr=abs(x-xO)xO=x;fvalue=subs(f,xO)dfvalue=subs(dfr xO);if (err<eps) | | (fvalue=0)/break,endend弦鼓法：secant, mfunction x,

35、iter,fvalue=secant(ff xOf xlf epszmaxiter)e弦截法X得到的近似解%iter迭代次数%fvalue函数在x处的值%f被求的非线性方程%x0, xl初始值%eps允许误差限%maxiter最大迭代次数fvalueO=subs(f,xO);fvalue=subs(f,xl);for iter=l:maxiterx=xl-fvalue*(xl-xO)/(fvalue-fvalue0)err=abs(x-xl)x0=xl;xl=x;fvalueO=subs(fz xO);fvalue=subs(fr xl)if (err<eps) | | (fvalue=

36、0)zbreakzendend求方程的实根：test8.mclear all：clc;syms xf=x.A3+2*x.A2+10*x-100;df=diff(fzxzl);eps=10e-6;maxiter=100;disp ( 1 Newton迭代初始值);xnl_0=0disp (1 Newton 迭代结果 1);xnlr iter_nlf fxnl=Newtoniter(f,df,xnl_0r epsrmaxiter)disp (* Newton迭代初始值1);xn2_0=5disp (1 Newton 迭代结果 1);xn2r iter_n2,fxn2=Newtoniter(f,df

37、,xn2_0r epsfmaxiter)disp.弦截法初始值,);卑微如螃蚁、坚强似大象xkl_O=Oxkl_l=ldisp (,弦截法迭代结果,)；xklz iter_klf fxkl=secant(fz xkl_Or xkl_l/epsrmaxiter)disp (,弦截法初始值,);xk2_0=5xk2_l=6disp弦截法迭代结果“xk2z iter_k2f fxk2=secant(f,xk2 0,xk2 1,epsrmaxiter)运行结果：Newton法结果:取两个不同初值0, 5kx(k)f|x(k)-x (k-1) |x(k)fI x (k) -x (k-1)i o * &#

38、187;0TOO- - - - -5125 * J;12 !101200103. 809522. 40581.1905!*i 2 *6. 5714335.86013. 42863. 48371.39060. 3258* 1 34. 546280. 75692. 02523. 46070. 00660. 0230! 43. 650811.82090. 89543. 46061.1043e-0041.5098e-007 ；i 53. 46770. 42770. 18303.4606-2. 8422e-0142. 5261e-009 jj 6 i3. 46066.3111e-0040. 007173

39、. 46061.3805e-0091.0559e-005:8 *1 »3. 4606-2.8422e-0142. 3098e-011 * i * 1 * i k*B*B*x(k)f|x(k)-x(k-1)|x(k)f|x(k) -x (k-1); *:0 .*0TOO 5125 * *9i 19 i .1-87162481.1905* * 27. 6923550. 43246. 69233. 983734. 80042.0163j 3B*IB1.9134-66. 53875. 77893. 654612.07110. 3291*| 4 99 92. 5366-45. 44240. 6

40、2323. 47981.15540. 1748 * : 内一3«I 5 .1B3. 879127. 25841. 34253. 46130. 04500. 0185 * * 63. 3758-4. 98050. 50343.46061.7917e-0047. 5046e-004| 73. 4535-0.42060. 07783.46062. 7963e-0082.9972e-006 j j 83. 45350. 00750. 0072：93. 4606-1.0939e-0051.2544e-004103. 4606-2.8388e-0101.8302e-0079.用二分法求方程/co

41、se + 2 = 0在区间欣4柯上的所有根.皿二分法,reseon. mfunction x,iter=resecm(f,a,b,eps)学二分法x近似解%iter迭代次数%f求解的方程%a, b求解区间%eps允许误差限fa=subs(f,a);fb=subs(f, b);iter=0;if(fa=O)x=a;returnendif (fb=O)x=b;returnendwhile(abs(a-b)>=eps)mf=subs(f,(a+b)/2);if(mf=O)x=mf;n=n+l;returnendif(mf*fa<0)b=(a+b)/2;elsea=(a+b)/2;endi

42、ter=iter+l;endx= (a+b)/2;iter=iter+l;求方程的实根：test9. mclear all;clc;syms xf=exp(x).*cos(x)+2;a=0;al=pi;a2=2*pi;a3=3*pi;b=4*pi;eps=10e-6;xl,iterl=resecm(f,a,alz eps)x2,iter2=resecm(fal,a2,eps)x3,iter3=resecm(fz &2,a3,eps)x4,iter4=resecm(fa3,b,eps)运行结果：O,pi区间的根x1 =1.8807； 迭代次数iterl = 20pi, 2*pi区间的根

43、x2 =4. 6941 ； 迭代次数 iter2 =202*pi,3*pi区间的根x3 =7. 8548；迭代次数 iter3 =203*pi,4*pi区间的根x4 =10. 9955；迭代次数 iter4 =2010.考虑函数f(2) = sin(7rr)? x 0,1,用等距节点作/(ar)的NeWo口插值,画出插值 多项式以及为的图像,观察收敛性.程序：N«wton 插值：Newtominter, mfunction f=Newtominter(x,y,xO)*牛顿插值x插值节点为对应的函数值e函数返回Newton插值多项式在x_0点的值fsyms t;if(length(x)

44、 = length(y)n = length(x);c(1:n) = 0.0;elsedisprx和y的维数不相等！,);return;endf = y(i)；yl = 0;1=1;for(i=l:n-l)for(j=i+l:n)yi(j) = (y(j)-y(i) )/(x(j)-x(i)； endc(i) = yl(i+l);1 = 1* (t-x(i);f = f + c(i)*l;simplify(f);y =山if(i=n-l)if (nargin = 3)*如果3个参数那么给出插值点的插值分f = subs(fz 1rx0);else+如果2个参数那么宜接给出插值多项式f = co

45、llect (f) ;e籽插值多项式展开f = vpa(fz 6);endend end 里微如螃蚁、坚强似大象用等距节点做f (x)的Newton插值：testlO. mnl=5;n2=10;n3=15;x0=0:0.01:l;yO=sin(pi.*x0);xl=Linspace (0, l,nl) ; %等距节点/ 节点数 5yl=sin(pi.*xl);f01=Newtominter(xlf ylf xO);x2=linspace (0,1, n2) ; %等儿节点,行点数 10y2=sin(pi.*x2);f02 = Newtominter(x2 r y2 r x0);x3=Linsp

46、ace (0,1, n3) ; +等距节点,节点数15y3=sin(pi.*x3);f03= Newtominter(x3 y3r xO);plot (xO, y0 1 -r 1) +1 京图hold onplot(x0, f01, '-g')%5 个节点plot(x0,f02r ,-k,)%10 个节点ploz(x0, f03,个节点legend (1原图1 J 5个节点Newton插值多项式,J ID个节点Newton插值多项式 J 15个 节点Newton插值多项式,)运行结果：取不同的节点做牛顿插值.得到结果图像如下：可以看出原图与插值多项式的图像近似重合,说明插值效果

47、较好.11.对函数幻=匕,±-5,5,取不同的节点数叫用等距节点作年尹近,插值, 尼察Runge现象.皿Lagrange 插值：Lagrange, mfunction fr fO = Lagrange(xfyfxO)%Lagrange插值x为插值结点,y为对应的函数值,xO为要计算的点.%函数返回L_n (x)表达式f和L_n (xO)的值f 0.syms t;if (Length(x) = length(y)n = length(x);elsedispx和y的维数不相等!,);return;ende检错f = 0.0;for(i = l:n)1 = y(i);for (j = l:

48、i-l)1 = 1* (t-x (j ) ) / (x(i) -x (j );end;for (j = i+l:n)1 = 1* (t-x (j)/(x(i)-x(j); *计算 Lagrange 基函数 end;f = f + 1;*计算Lagrange插值函数simplify (f) ;+化简if (i=n)if(nargin = 3)f0 = subs(fr;*如果3个参数那么计算插值点的函数值elsef = collect (f) ;*如果2个参数那么将插值多项式展开f = vpa(f,6);e将插值多项式的系数化成6位精度的小数endendend用等距节点做Lagrange插值：te

49、st, mclear all;clc;nl=5;n2=10;n3=15;x0=-5:0.02:5;yO=L./(l+xO.A2);xl=linspace (-5, 5, nl) ; +等距节点,节点数 5yl=l./(1+xl.A2);f1,f01 = Lagrange(xlr ylr xO);x2=linspace (-5,5r n2) ; %等距节点/ 节点数 10y2=L./(l+x2.A2);f2,f02 = Lagrange(x2,y2,xO);x3=Linspace (-5, 5, n3) ; *等距节点/ 节点数 15y3=l./(l+x3.A2);f3,f03 = Lagrange(x3r y3z xO);plot (xO, y0 1 -r 1) +1 取图hold onplot(x0,f01f ,-b,)%5 个节点plot(xOr f02r l-gl)%10 个节点ploz(x0, f03r ,-kl)%15 个节点xlabel (1 x 1) ; ylabel (1 f (x) r L (x) 1);legend (1原图1 / 1 5个节点Lagrange插值多项式1 / 110个节点Lagrange插值多项式 ,15个节点Lagrange插值多项式,)运行结果：选取了 5. 10. 15个节点做Lagrange插值,得到原图与