大連理工大學矩陣與數(shù)值分析上機作業(yè)

1、矩陣與數(shù)值分析上機作業(yè) 學校： 大連理工大學 學院： 班級： 姓名： 學號： 授課老師： 注：編程語言Matlab程序：Norm.m函數(shù)function s=Norm(x,m)%求向量x的范數(shù)%m取1,2,inf分別 表示1,2，無窮范數(shù)n=length(x);s=0;switch m case 1 %1-范數(shù) for i=1:n s=s+abs(x(i); end case 2 %2-范數(shù) for i=1:n s=s+x(i)2; end s=sqrt(s); case inf %無窮-范數(shù) s=max(abs(x);end 計算向量x，y的范數(shù)Test1.mclear all;clc;n1

2、=10;n2=100;n3=1000;x1=1./1:n1'x2=1./1:n2'x3=1./1:n3'y1=1:n1'y2=1:n2'y3=1:n3'disp('n=10時');disp('x的1-范數(shù):');disp(Norm(x1,1);disp('x的2-范數(shù):');disp(Norm(x1,2);disp('x的無窮-范數(shù):');disp(Norm(x1,inf);disp('y的1-范數(shù):');disp(Norm(y1,1);disp('y的2-范

3、數(shù):');disp(Norm(y1,2);disp('y的無窮-范數(shù):');disp(Norm(y1,inf);disp('n=100時');disp('x的1-范數(shù):');disp(Norm(x2,1);disp('x的2-范數(shù):');disp(Norm(x2,2);disp('x的無窮-范數(shù):');disp(Norm(x2,inf);disp('y的1-范數(shù):');disp(Norm(y2,1);disp('y的2-范數(shù):');disp(Norm(y2,2);disp(&

4、#39;y的無窮-范數(shù):');disp(Norm(y2,inf);disp('n=1000時');disp('x的1-范數(shù):');disp(Norm(x3,1);disp('x的2-范數(shù):');disp(Norm(x3,2);disp('x的無窮-范數(shù):');disp(Norm(x3,inf);disp('y的1-范數(shù):');disp(Norm(y3,1);disp('y的2-范數(shù):');disp(Norm(y3,2);disp('y的無窮-范數(shù):');disp(Norm(y

5、3,inf);運行結果：n=10時x的1-范數(shù):2.9290；x的2-范數(shù):1.2449； x的無窮-范數(shù):1y的1-范數(shù):55； y的2-范數(shù):19.6214； y的無窮-范數(shù):10n=100時x的1-范數(shù):5.1874；x的2-范數(shù): 1.2787； x的無窮-范數(shù):1y的1-范數(shù):5050； y的2-范數(shù):581.6786； y的無窮-范數(shù):100n=1000時x的1-范數(shù):7.4855； x的2-范數(shù):1.2822； x的無窮-范數(shù):1y的1-范數(shù): 500500； y的2-范數(shù):1.8271e+004；y的無窮-范數(shù):1000程序Test2.mclear all;clc;n=100;

6、%區(qū)間h=2*10(-15)/n;%步長x=-10(-15):h:10(-15);%第一種原函數(shù)f1=zeros(1,n+1);for k=1:n+1 if x(k)=0 f1(k)=log(1+x(k)/x(k); else f1(k)=1; endendsubplot(2,1,1);plot(x,f1,'-r');axis(-10(-15),10(-15),-1,2);legend('原圖');%第二種算法f2=zeros(1,n+1);for k=1:n+1 d=1+x(k); if(d=1) f2(k)=log(d)/(d-1); else f2(k)=

7、1; endendsubplot(2,1,2);plot(x,f2,'-r');axis(-10(-15),10(-15),-1,2);legend('第二種算法');運行結果：顯然第二種算法結果不準確，是因為計算機中的舍入誤差造成的，當時，計算機進行舍入造成恒等于1，結果函數(shù)值恒為1。程序：秦九韶算法:QinJS.mfunction y=QinJS(a,x)%y輸出函數(shù)值%a多項式系數(shù)，由高次到零次%x給定點n=length(a);s=a(1);for i=2:n s=s*x+a(i);endy=s;計算p(x):test3.mclear all;clc;x=

8、1.6:0.2:2.4;%x=2的鄰域disp('x=2的鄰域：');xa=1 -18 144 -672 2016 -4032 5376 -4608 2304 -512;p=zeros(1,5);for i=1:5 p(i)=QinJS(a,x(i);enddisp('相應多項式p值：');pxk=1.95:0.01:20.5;nk=length(xk);pk=zeros(1,nk);k=1;for k=1:nk pk(k)=QinJS(a,xk(k);endplot(xk,pk,'-r');xlabel('x');ylabel(

9、'p(x)');運行結果：x=2的鄰域：x =1.6000 1.8000 2.0000 2.2000 2.4000相應多項式p值：p = 1.0e-003 * -0.2621 -0.0005 0 0.0005 0.2621p(x)在1.95,20.5上的圖像程序：LU分解，LUDe.mfunction L,U=LUDe.(A)%不帶列主元的LU分解N = size(A);n = N(1);L=eye(n);U=zeros(n);for i=1:n U(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);endfor i=2:n for j=i:n z=0; for

10、 k=1:i-1 z=z+L(i,k)*U(k,j); end U(i,j)=A(i,j)-z; end for j=i+1:n z=0; for k=1:i-1 z=z+L(j,k)*U(k,i); end L(j,i)=(A(j,i)-z)/U(i,i); endendPLU分解，PLUDe.mfunction P,L,U =PLUDe.(A)%帶列主元的LU分解m,m=size(A);U=A;P=eye(m);L=eye(m);for i=1:m for j=i:m t(j)=U(j,i); for k=1:i-1 t(j)=t(j)-U(j,k)*U(k,i); end end a=i

11、;b=abs(t(i); for j=i+1:m if b<abs(t(j) b=abs(t(j); a=j; end end if a=i for j=1:m c=U(i,j); U(i,j)=U(a,j); U(a,j)=c; end for j=1:m c=P(i,j); P(i,j)=P(a,j); P(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end U(i,i)=t(i); for j=i+1:m U(j,i)=t(j)/t(i); end for j=i+1:m for k=1:i-1 U(i,j)=U(i,j)-U(i,k)*U(k,j

12、); end endendL=tril(U,-1)+eye(m);U=triu(U,0);（1） （2）程序：Test4.mclear all;clc;for n=5:30 x=zeros(n,1); A=-ones(n); A(:,n)=ones(n,1); for i=1:n A(i,i)=1; for j=(i+1):(n-1) A(i,j)=0; end x(i)=1/i; end disp('當n=');disp(n); disp('方程精確解:'); x b=A*x; %系數(shù)b disp('利用LU分解方程組的解:'); L,U=LU

13、De.(A); %LU分解 xLU=U(Lb) disp('利用PLU分解方程組的解:'); P,L,U =PLUDe.(A); %PLU分解 xPLU=U(L(Pb) %求解A的逆矩陣 disp('A的準確逆矩陣:'); InvA=inv(A) InvAL=zeros(n); %利用LU分解求A的逆矩陣 I=eye(n); for i=1:n InvAL(:,i)=U(LI(:,i); end disp('利用LU分解的A的逆矩陣:'); InvALEnd運行結果：（1） 只列出n=5,6,7的結果當n= 5方程精確解:x =1.0000 0.

14、5000 0.3333 0.2500 0.2000利用LU分解方程組的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000利用PLU分解方程組的解:xPLU = 1.0000 0.5000 0.3333 0.2500 0.2000當n=6方程精確解:x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667利用LU分解方程組的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667利用PLU分解方程組的解:xPLU = 1.0000 0.5000 0.3333 0.2500 0.2000

15、0.1667當n= 7方程精確解:x = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429利用LU分解方程組的解:xLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429利用PLU分解方程組的解:xPLU = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429（2） 只列出n=5,6,7時A的逆矩陣的結果當n= 5A的準確逆矩陣:InvA = 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.5000 -0.250

16、0 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0625利用LU分解的A的逆矩陣:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0.5000 -0.2500 -0.2500 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625當n= 6A的準確逆矩陣:InvA = 0.5000

17、 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313利用LU分解的A的逆矩陣:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.12

18、50 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0313當n= 7A的準確逆矩陣:InvA = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0.1250 -0.

19、0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156利用LU分解的A的逆矩陣:InvAL = 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0 0.5000 -0.2500 -0

20、.1250 -0.0625 -0.0625 0 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156程序：Cholesky分解：Cholesky.mfunction L=Cholesky(A)N = size(A);n = N(1);L=zeros(n);L(1,1)=sqrt(A(1,1);for i=2:n L(i,1)=A(i,1)/L(1,1);endfor j=2

21、:n s1=0; for k=1:j-1 s1=s1+L(j,k)2; end L(j,j)=sqrt(A(j,j)-s1); for i=j+1:n s2=0; for k=1:j-1 s2=s2+L(i,k)*L(j,k); end L(i,j)=(A(i,j)-s2)/L(j,j); endend計算Ax=b；Test5.mclear all;clc;for n=10:20 A=zeros(n,n); b=zeros(n,1); for i=1:n for j=1:n A(i,j)=1/(i+j-1); end b(i,1)=i; end disp('n=');disp(

22、n); disp('方程組原始解');x0=Ab disp('利用Cholesky分解的方程組的解'); L=Cholesky(A) x=L'(Lb)end運行結果：只列出了n=10,11的結果n=10方程組原始解x0 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2108 3.5947 -6.3233 6.5114 -3.6233 0.8407利用Cholesky分解的方程組的解x = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2105 3.5939 -6.32

23、19 6.5100 -3.6225 0.8405n= 11方程組原始解x0 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4039 3.2863 -4.7869 4.2260 -2.0685 0.4305利用Cholesky分解的方程組的解x = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0563 0.3668 -1.3972 3.2716 -4.7669 4.2094 -2.06080.4290程序：（1）House.mfunction u=House(x)n=length(x);e1=eye(n,1)

24、;w=x-norm(x,2)*e1;u=w/norm(w,2);(2)Hou_A.mfunction HA=Hou_A(A)a1=A(:,1);n=length(a1);e1=eye(n,1);w=a1-norm(a1,2)*e1;u=w/norm(w,2);H=eye(n)-2*u*u'HA=H*A;(3)test6.mclear all;clc;A=1 2 3 4; -1 2 sqrt(2) sqrt(3); -2 2 exp(1) pi; -sqrt(10) 2 -3 7; 0 2 7 5/2;HA=Hou_A(A)運行結果：H = 0.2500 -0.2500 -0.5000

25、-0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0.5270 0.1667 0 0 0 0 0 1.0000HA = 4.0000 -2.5811 1.4090 -6.5378 0.0000 0.4730 0.8839 -1.7805 0.0000 -1.0541 1.6576 -3.8836 0.0000 -2.8289 -4.6770 -4.1078 0 2.0000 7.0000 2.5000程序：Jacobi迭代：Jaccobi.mfunction

26、x,n=Jaccobi(A,b,x0)%-·方程組系數(shù)陣A%-·方程組右端頂b%- 初始值x0%-求解要求精確度eps%-迭代步數(shù)控制M%-·返回求得的解x %-·返回迭代步數(shù)nM=1000;eps=1.0e-5;D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣J=D(L+U);f=Db;x=J*x0+fn=1; %迭代次數(shù)err=norm(x-x0,inf)while(err>=eps) x0=x; x=J*x0+f n=n+1; err=norm(x-x0

27、,inf) if(n>=M) disp('Warning: 迭代次數(shù)太多，可能不收斂?'); return; endendGauss_Seidel迭代：Gauss_Seidel.mfunction x,n=Gauss_Seidel(A,b,x0)%-Gauss-Seidel迭代法解線性方程組%-方程組系數(shù)陣 A%-方程組右端項 b%-初始值 x0%-求解要求的精確度 eps%-迭代步數(shù)控制 M%-返回求得的解 x %-返回迭代步數(shù) neps=1.0e-5;M=10000;D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-t

28、riu(A,1); %求A的上三角陣G=(D-L)U;f=(D-L)b;x=G*x0+fn=1; %迭代次數(shù)err=norm(x-x0,inf)while(err>=eps) x0=x; x=G*x0+f n=n+1; err=norm(x-x0,inf) if(n>=M) disp('Warning: 迭代次數(shù)太多，可能不收斂！'); return; endend解方程組，test7.mclear all;clc;A=5 -1 -3; -1 2 4; -3 4 15;b=-2;1;10;disp('精確解');x=Abdisp('迭代初始值

29、');x0=0;0;0disp('Jacobi迭代過程：');xj,nj=Jaccobi(A,b,x0);disp('Jacobi最終迭代結果：');xjdisp('迭代次數(shù)');njdisp('Gauss-Seidel迭代過程：');xg,ng=Gauss_Seidel(A,b,x0);disp('Gauss-Seidel最終迭代結果：');xgdisp('迭代次數(shù)');ng運行結果：精確解x = -0.0820 -1.8033 1.1311迭代初始值x0 = 0 0 0Jacobi迭代過

30、程：x = -0.4000 0.5000 0.6667err = 0.6667x = 0.1000 -1.0333 0.4533err =1.5333.x = -0.0820 -1.8033 1.1311err = 9.6603e-006Jacobi最終迭代結果：xj = -0.0820 -1.8033 1.1311迭代次數(shù)nj = 281Gauss-Seidel迭代過程：x = -0.4000 0.3000 0.5067err = 0.5067x = -0.0360 -0.5313 0.8012err = 0.8313x = -0.0256 -1.1151 0.9589err =0.5838

31、.x = -0.0820 -1.8033 1.1311err = 9.4021e-006Gauss-Seidel最終迭代結果：xg = -0.0820 -1.8033 1.1311迭代次數(shù)ng =20程序：Newton迭代法：Newtoniter.mfunction x,iter,fvalue=Newtoniter(f,df,x0,eps,maxiter)%牛頓法 x得到的近似解%iter迭代次數(shù)%fvalue函數(shù)在x處的值%f，df被求的非線性方程及導函數(shù)%x0初始值%eps 允許誤差限%maxiter 最大迭代次數(shù)fvalue=subs(f,x0);dfvalue=subs(df,x0);

32、for iter=1:maxiter x=x0-fvalue/dfvalue err=abs(x-x0) x0=x; fvalue=subs(f,x0) dfvalue=subs(df,x0); if(err<eps)|(fvalue=0),break,endend弦截法：secant.mfunction x,iter,fvalue=secant(f,x0,x1,eps,maxiter)%弦截法 x得到的近似解%iter迭代次數(shù)%fvalue函數(shù)在x處的值%f被求的非線性方程%x0，x1初始值%eps 允許誤差限%maxiter 最大迭代次數(shù)fvalue0=subs(f,x0);fval

33、ue=subs(f,x1);for iter=1:maxiter x=x1-fvalue*(x1-x0)/(fvalue-fvalue0) err=abs(x-x1) x0=x1;x1=x; fvalue0=subs(f,x0);fvalue=subs(f,x1) if(err<eps)|(fvalue=0),break,endend求方程的實根：test8.mclear all；clc;syms xf=x.3+2*x.2+10*x-100;df=diff(f,x,1);eps=10e-6;maxiter=100;disp('Newton迭代初始值');xn1_0=0di

34、sp('Newton迭代結果');xn1,iter_n1,fxn1=Newtoniter(f,df,xn1_0,eps,maxiter)disp('Newton迭代初始值');xn2_0=5disp('Newton迭代結果');xn2,iter_n2,fxn2=Newtoniter(f,df,xn2_0,eps,maxiter)disp('弦截法初始值');xk1_0=0xk1_1=1disp('弦截法迭代結果');xk1,iter_k1,fxk1=secant(f,xk1_0,xk1_1,eps,maxiter)

35、disp('弦截法初始值');xk2_0=5 xk2_1=6disp('弦截法迭代結果');xk2,iter_k2,fxk2=secant(f,xk2_0,xk2_1,eps,maxiter) 運行結果：Newton法結果:取兩個不同初值0，5kx(k)f|x(k)-x(k-1)|x(k)f|x(k)-x(k-1)00-10051251101200103.809522.40581.190526.5714335.86013.42863.48371.39060.325834.546280.75692.02523.46070.00660.023043.650811.8

36、2090.89543.46061.1043e-0041.5098e-00753.46770.42770.18303.4606-2.8422e-0142.5261e-00963.46066.3111e-0040.007173.46061.3805e-0091.0559e-00583.4606-2.8422e-0142.3098e-011弦截法迭代結果：取兩種不同初值0,1；5,6kx(k)f|x(k)-x(k-1)|x(k)f|x(k)-x(k-1)00-100512511-87162481.190527.6923550.43246.69233.983734.80042.016331.9134-

37、66.5387 5.77893.654612.07110.329142.5366-45.44240.62323.47981.15540.174853.879127.25841.34253.46130.04500.018563.3758-4.98050.50343.46061.7917e-0047.5046e-00473.4535-0.42060.07783.46062.7963e-0082.9972e-00683.45350.00750.007293.4606-1.0939e-0051.2544e-004103.4606-2.8388e-0101.8302e-007程序：二分法：resecm.

38、mfunction x,iter=resecm(f,a,b,eps)%二分法 x 近似解%iter 迭代次數(shù)%f 求解的方程%a，b 求解區(qū)間%eps 允許誤差限fa=subs(f,a);fb=subs(f,b);iter=0;if(fa=0) x=a;returnendif(fb=0) x=b;returnendwhile(abs(a-b)>=eps) mf=subs(f,(a+b)/2); if(mf=0) x=mf;n=n+1;return end if(mf*fa<0) b=(a+b)/2; else a=(a+b)/2; end iter=iter+1;endx=(a+b

39、)/2;iter=iter+1;求方程的實根：test9.mclear all;clc;syms xf=exp(x).*cos(x)+2;a=0;a1=pi;a2=2*pi;a3=3*pi;b=4*pi;eps=10e-6;x1,iter1=resecm(f,a,a1,eps)x2,iter2=resecm(f,a1,a2,eps)x3,iter3=resecm(f,a2,a3,eps)x4,iter4=resecm(f,a3,b,eps) 運行結果：0,pi區(qū)間的根x1 =1.8807； 迭代次數(shù)iter1 = 20pi,2*pi區(qū)間的根x2 =4.6941； 迭代次數(shù)iter2 =202*

40、pi,3*pi區(qū)間的根x3 =7.8548； 迭代次數(shù)iter3 =203*pi,4*pi區(qū)間的根x4 =10.9955；迭代次數(shù)iter4 =20程序：Newton插值：Newtominter.mfunction f=Newtominter(x,y,x0)%牛頓插值 x插值節(jié)點%y為對應的函數(shù)值%函數(shù)返回Newton插值多項式在x_0點的值fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp('x和y的維數(shù)不相等！'); return;end f = y(1);y1 = 0;l = 1;f

41、or(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f); y = y1; if(i=n-1) if(nargin = 3) %如果3個參數(shù)則給出插值點的插值結果 f = subs(f,'t',x0); else %如果2個參數(shù)則直接給出插值多項式 f = collect(f); %將插值多項式展開 f = vpa(f, 6); end endend用等距節(jié)點做f(x)的Newton插值：test10

42、.mn1=5;n2=10;n3=15;x0=0:0.01:1;y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距節(jié)點,節(jié)點數(shù)5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace(0,1,n2);%等距節(jié)點,節(jié)點數(shù)10y2=sin(pi.*x2);f02 = Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距節(jié)點,節(jié)點數(shù)15y3=sin(pi.*x3);f03= Newtominter(x3,y3,x0);plot(x0,y0,'-r')%原圖hold onpl

43、ot(x0,f01,'-g')%5個節(jié)點plot(x0,f02,'-k')%10個節(jié)點plot(x0,f03,'-b')%15個節(jié)點legend('原圖','5個節(jié)點Newton插值多項式','10個節(jié)點Newton插值多項式','15個節(jié)點Newton插值多項式')運行結果：取不同的節(jié)點做牛頓插值。得到結果圖像如下：可以看出原圖與插值多項式的圖像近似重合，說明插值效果較好。程序：Lagrange插值：Lagrange.mfunction f,f0 = Lagrange(x,y,x0

44、) %Lagrange插值 x為插值結點，y為對應的函數(shù)值，x0為要計算的點。%函數(shù)返回L_n(x)表達式f和L_n(x0)的值f0。syms t;if(length(x) = length(y) n = length(x); else disp('x和y的維數(shù)不相等！'); return;end %檢錯 f = 0.0;for(i = 1:n) l = y(i); for(j = 1:i-1) l = l*(t-x(j)/(x(i)-x(j); end; for(j = i+1:n) l = l*(t-x(j)/(x(i)-x(j); %計算Lagrange基函數(shù) end;

45、f = f + l; %計算Lagrange插值函數(shù) simplify(f); %化簡 if(i=n) if(nargin = 3) f0 = subs(f,'t',x0); %如果3個參數(shù)則計算插值點的函數(shù)值 else f = collect(f); %如果2個參數(shù)則將插值多項式展開 f = vpa(f,6); %將插值多項式的系數(shù)化成6位精度的小數(shù) end endend用等距節(jié)點做Lagrange插值：test11.mclear all;clc;n1=5;n2=10;n3=15;x0=-5:0.02:5;y0=1./(1+x0.2);x1=linspace(-5,5,n1);%等距節(jié)點,節(jié)點數(shù)5y1=1./(1+x1.2);f1,f01 = Lagrange(x1,y1,x0);x2=linspace(-5,5,n2);%等距節(jié)點,節(jié)點數(shù)10y2=1./(1+x2.2);f2,f02 = Lagrange(x2,y2,x0);x3=linspace(-5,5,n3);%等距節(jié)點,節(jié)點數(shù)15y3=1./(1+x3.2);f3,f03 = Lagrange(x3,y3,x0);plot(x0,y0,'-r')