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

矩陣與數(shù)值分析上機(jī)作業(yè)學(xué)院：班級(jí)：姓名：學(xué)號(hào)：授課老師：注：編程語(yǔ)言Matlab程序：Norm.m函數(shù)functions=Norm(x,m)%x的范數(shù)%m1,2,inf1,2，無(wú)窮范數(shù)n=length(x);s=0;switchmcase1 %1-fori=1:ns=s+abs(x(i));endcase2 %2-fori=1:ns=s+x(i)^2;ends=sqrt(s);caseinf%無(wú)窮-范數(shù)s=max(abs(x));endx，yTest1.mclearall;clc;n1=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時(shí)');disp('x1-范數(shù):');disp(Norm(x1,1));disp('x2-范數(shù):');disp(Norm(x1,2));disp('x的無(wú)窮-范數(shù):');disp(Norm(x1,inf));disp('y1-范數(shù):');disp(Norm(y1,1));disp('y2-范數(shù):');disp(Norm(y1,2));disp('y的無(wú)窮-范數(shù):');disp(Norm(y1,inf));disp('n=100時(shí)');disp('x1-范數(shù):');disp(Norm(x2,1));disp('x2-范數(shù):');disp(Norm(x2,2));disp('x的無(wú)窮-范數(shù):');disp(Norm(x2,inf));disp('y1-范數(shù):');disp(Norm(y2,1));disp('y2-范數(shù):');disp(Norm(y2,2));disp('y的無(wú)窮-范數(shù):');disp(Norm(y2,inf));disp('n=1000時(shí)');disp('x1-范數(shù):');disp(Norm(x3,1));disp('x2-范數(shù):');disp(Norm(x3,2));disp('x的無(wú)窮-范數(shù):');disp(Norm(x3,inf));disp('y1-范數(shù):');disp(Norm(y3,1));disp('y2-范數(shù):');disp(Norm(y3,2));disp('y的無(wú)窮-范數(shù):');disp(Norm(y3,inf));運(yùn)行結(jié)果：n=10時(shí)x的1-范:2.9290；x的2-范:1.2449；x的無(wú)-范:1y的1-范:55； y的2-范:19.6214；y的無(wú)-范n=100時(shí)x1:5.1874；x2-范數(shù):1.2787；x:1y1:5050；y2-:581.6786；y:100n=1000時(shí)x的1-范:7.4855；x的2-范數(shù):1.2822； x的無(wú)范:1y1500500；y2-范數(shù):1.8271e+004；y:1000程序Test2.mclearclc;h=2*10^(-15)/n;%步長(zhǎng)x=-10^(-15):h:10^(-15);%第一種原函數(shù)f1=zeros(1,n+1);fork=1:n+1ifx(k)~=0f1(k)=log(1+x(k))/x(k);elsef1(k)=1;endendsubplot(2,1,1);plot(x,f1,'-r');axis([-10^(-15),10^(-15),-1,2]);legend('原圖');%第二種算法f2=zeros(1,n+1);fork=1:n+1d=1+x(k);if(d~=1)f2(k)=log(d)/(d-1);elsef2(k)=1;endendsubplot(2,1,2);plot(x,f2,'-r');axis([-10^(-15),10^(-15),-1,2]);legend('第二種算法');運(yùn)行結(jié)果：顯然第二種算法結(jié)果不準(zhǔn)確計(jì)算機(jī)進(jìn)行舍入造成d11。程序：秦九韶算法:QinJS.mfunctiony=QinJS(a,x)%y輸出函數(shù)值%a多項(xiàng)式系數(shù)，由高次到零次%x給定點(diǎn)n=length(a);s=a(1);fori=2:ns=s*x+a(i);endy=s;計(jì)算p(x):test3.mclearall;clc;x=1.6:0.2:2.4;%x=2的鄰域disp('x=2的鄰域：');xa=[1-18144-6722016-40325376-46082304-512];p=zeros(1,5);fori=1:5p(i)=QinJS(a,x(i));enddisp('p值：');pxk=1.95:0.01:20.5;nk=length(xk);pk=zeros(1,nk);k=1;fork=1:nkpk(k)=QinJS(a,xk(k));endplot(xk,pk,'-r');xlabel('x');ylabel('p(x)');運(yùn)行結(jié)果：x=2的鄰域：x=1.6000 1.8000 2.0000 2.2000 2.4000相應(yīng)多項(xiàng)式p值：p=1.0e-003*-0.2621 -0.0005 0 0.0005 0.2621p(x)x[1.95,20.5]上的圖像程序：LU分解，LUDecom.mfunction[L,U]=LUDecom(A)%LU分解N=size(A);n=N(1);fori=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);endfori=2:nforj=i:nz=0;fork=1:i-1z=z+L(i,k)*U(k,j);endU(i,j)=A(i,j)-z;endforj=i+1:nz=0;fork=1:i-1z=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分解[m,m]=size(A);U=A;P=eye(m);L=eye(m);fori=1:mforj=i:mfork=1:i-1t(j)=t(j)-U(j,k)*U(k,i);endenda=i;b=abs(t(i));forj=i+1:mifb<abs(t(j))b=abs(t(j));a=j;endendifa~=iforj=1:mc=U(i,j);U(i,j)=U(a,j);U(a,j)=c;endforj=1:mc=P(i,j);P(i,j)=P(a,j);P(a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endU(i,i)=t(i);forj=i+1:mU(j,i)=t(j)/t(i);endforj=i+1:mfork=1:i-1U(i,j)=U(i,j)-U(i,k)*U(k,j);endendendL=tril(U,-1)+eye(m);U=triu(U,0);（1）（2）程序：Test4.mclearall;clc;forn=5:30x=zeros(n,1);A=-ones(n);A(:,n)=ones(n,1);fori=1:nA(i,i)=1;forj=(i+1):(n-1)A(i,j)=0;endx(i)=1/i;enddisp('當(dāng)n=');disp(n);disp('方程精確解:');xb=A*x; %bLU分解方程組的解:');[L,U]=LUDecom(A);%LU分解xLU=U\(L\b)利用PLU分解方程組的解:');[P,L,U]=PLUDecom(A); %PLU分xPLU=U\(L\(P\b))%A的逆矩陣的準(zhǔn)確逆矩陣:');InvA=inv(A)InvAL=zeros(n)LUA的逆矩陣I=eye(n);fori=1:nInvAL(:,i)=U\(L\I(:,i));endLUA的逆矩陣:');InvALEnd運(yùn)行結(jié)果：n=5,6,7當(dāng)n=5方程精確解:x=1.00000.50000.33330.25000.2000LUxLU=1.00000.50000.33330.25000.2000PLUxPLU=1.00000.50000.33330.25000.2000當(dāng)n=6方程精確解:x=1.00000.50000.33330.25000.20000.1667LUxLU=1.00000.50000.33330.25000.20000.1667PLUxPLU=1.00000.50000.33330.25000.20000.1667當(dāng)n=7方程精確解:x=1.00000.50000.33330.25000.20000.16670.1429LUxLU=1.00000.50000.33330.25000.20000.16670.1429PLUxPLU=1.00000.50000.33330.25000.20000.16670.1429n=5,6,7A當(dāng)n=5A的準(zhǔn)確逆矩陣:InvA=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.0625LU分解的AInvAL=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當(dāng)n=6A的準(zhǔn)確逆矩陣:InvA=0.5000-0.2500-0.1250-0.0625-0.0313-0.031300.5000-0.2500-0.1250-0.0625-0.0625000.5000-0.2500-0.1250-0.12500000.5000-0.2500-0.250000000.5000-0.50000.50000.25000.12500.06250.03130.0313LU分解的AInvAL=0.5000-0.2500-0.1250-0.0625-0.0313-0.031300.5000-0.2500-0.1250-0.0625-0.0625000.5000-0.2500-0.1250-0.12500000.5000-0.2500-0.250000000.5000-0.50000.50000.25000.12500.06250.03130.0313當(dāng)n=7A的準(zhǔn)確逆矩陣: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.01560.0156LU分解的AInvAL=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.01560.0156程序：Cholesky分解：Cholesky.mfunctionL=Cholesky(A)N=size(A);n=N(1);L=zeros(n);L(1,1)=sqrt(A(1,1));fori=2:nL(i,1)=A(i,1)/L(1,1);endforj=2:ns1=0;fork=1:j-1s1=s1+L(j,k)^2;endL(j,j)=sqrt(A(j,j)-s1);fori=j+1:ns2=0;fork=1:j-1s2=s2+L(i,k)*L(j,k);endL(i,j)=(A(i,j)-s2)/L(j,j);endend計(jì)算Ax=b；Test5.mclearall;clc;forn=10:20A=zeros(n,n);fori=1:nforj=1:nA(i,j)=1/(i+j-1);endb(i,1)=i;end方程組原始解');x0=A\bCholesky分解的方程組的解');L=Cholesky(A)x=L'\(L\b)end運(yùn)行結(jié)果：只列出了n=10,11n=10x0=1.0e+008*-0.00000.0010-0.02330.2330-1.21083.5947-6.32336.5114-3.62330.8407利用Cholesky分解的方程組的解x=1.0e+008*-0.00000.0010-0.02330.2330-1.21053.5939-6.32196.5100-3.62250.8405n=11方程組原始解x0=1.0e+009*0.0000-0.00020.0046-0.05670.3687-1.40393.2863-4.78694.2260-2.06850.4305利用Cholesky分解的方程組的解x=1.0e+009*0.0000-0.00020.0046-0.05630.3668-1.39723.2716-4.76694.2094-2.06080.4290程序：（1）House.mfunctionu=House(x)n=length(x);e1=eye(n,1);w=x-norm(x,2)*e1;u=w/norm(w,2);(2)Hou_A.mfunctionHA=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.mclearall;clc;A=[1234;-12sqrt(2)sqrt(3);-22exp(1)pi;-sqrt(10)2-37;0275/2];HA=Hou_A(A)運(yùn)行結(jié)果：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.0000HA=4.0000-2.58111.4090-6.53780.00000.47300.8839-1.78050.0000-1.05411.6576-3.88360.0000-2.8289-4.6770-4.107802.00007.00002.5000程序：Jacobi迭代：Jaccobi.mfunction[x,n]=Jaccobi(A,b,x0)%--·A%--·b%--初始值x0%--求解要求精確度eps%--迭代步數(shù)控制M%--·返回求得的解x%--·返回迭代步數(shù)nM=1000;eps=1.0e-5;D=diag(diag(A));A的對(duì)角矩陣L=-tril(A,-1);A的下三角陣U=-triu(A,1);A的上三角陣J=D\(L+U);f=D\b;x=J*x0+fn=1;%迭代次數(shù)err=norm(x-x0,inf)while(err>=eps)x0=x;x=J*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M):迭代次數(shù)太多，可能不收斂?'return;endendGauss_Seidel迭代：Gauss_Seidel.mfunction[x,n]=Gauss_Seidel(A,b,x0)%--Gauss-Seidel迭代法解線性方程組%--方程組系數(shù)陣A%--方程組右端項(xiàng)b%--初始值x0%--求解要求的精確度eps%--迭代步數(shù)控制M%--返回求得的解x%--返回迭代步數(shù)neps=1.0e-5;M=10000;D=diag(diag(A));A的對(duì)角矩陣L=-tril(A,-1);A的下三角陣U=-triu(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+fn=n+1;err=norm(x-x0,inf)if(n>=M):迭代次數(shù)太多，可能不收斂！'return;endend解方程組，test7.mclearall;clc;A=[5-1-3;-124;-3415];b=[-2;1;10];disp('精確解');x=A\bdisp('迭代初始值');x0=[0;0;0]disp('Jacobi迭代過(guò)程：');[xj,nj]=Jaccobi(A,b,x0);disp('Jacobi最終迭代結(jié)果：');xjdisp('迭代次數(shù)');njdisp('Gauss-Seidel迭代過(guò)程：');[xg,ng]=Gauss_Seidel(A,b,x0);disp('Gauss-Seidel最終迭代結(jié)果：');xgdisp('迭代次數(shù)');ng運(yùn)行結(jié)果：精確解x=-0.0820-1.80331.1311迭代初始值x0=000Jacobi迭代過(guò)程：x=-0.40000.50000.6667err=0.6667x=0.1000-1.03330.4533err=1.5333...x=-0.0820-1.80331.1311err=9.6603e-006Jacobi最終迭代結(jié)果：xj=-0.0820-1.80331.1311迭代次數(shù)nj=281Gauss-Seidel迭代過(guò)程：x=-0.40000.30000.5067err=0.5067x=-0.0360-0.53130.8012err=0.8313x=-0.0256-1.11510.9589err=0.5838...x=-0.0820-1.80331.1311err=9.4021e-006Gauss-Seidel最終迭代結(jié)果：xg=-0.0820-1.80331.1311迭代次數(shù)ng=20程序：Newton迭代法：Newtoniter.mfunction[x,iter,fvalue]=Newtoniter(f,df,x0,eps,maxiter)%牛頓法x得到的近似解%iter迭代次數(shù)%fvaluex處的值%f，df被求的非線性方程及導(dǎo)函數(shù)%x0初始值%eps允許誤差限%maxiter最大迭代次數(shù)fvalue=subs(f,x0);dfvalue=subs(df,x0);foriter=1:maxiterx=x0-fvalue/dfvalueerr=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ù)%fvaluex處的值%f被求的非線性方程%x0，x1初始值%eps允許誤差限%maxiter最大迭代次數(shù)fvalue0=subs(f,x0);fvalue=subs(f,x1);foriter=1:maxiterx=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求方程的實(shí)根：test8.mclearall；clc;symsxf=x.^3+2*x.^2+10*x-100;df=diff(f,x,1);eps=10e-6;maxiter=100;disp('Newton迭代初始值');xn1_0=0disp('Newton迭代結(jié)果');disp('Newton迭代初始值');xn2_0=5disp('Newton迭代結(jié)果');[xn2,iter_n2,fxn2]=Newtoniter(f,df,xn2_0,eps,maxiter)disp('弦截法初始值');xk1_0=0xk1_1=1disp('弦截法迭代結(jié)果');[xk1,iter_k1,fxk1]=secant(f,xk1_0,xk1_1,eps,maxiter)disp('弦截法初始值');xk2_0=5xk2_1=6disp('弦截法迭代結(jié)果');[xk2,iter_k2,fxk2]=secant(f,xk2_0,xk2_1,eps,maxiter)運(yùn)行結(jié)果：Newton法結(jié)果:取兩個(gè)不同初值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.82090.89543.46061.1043e-0041.5098e-00753.46770.42770.18303.4606-2.8422e-0142.5261e-00963.46066.3111e-0040.007173.46061.3805e-0091.0559e-005883.4606-2.8422e-0142.3098e-011弦截法迭代結(jié)果：取兩種不同初值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-66.53875.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.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)endif(mf*fa<0)b=(a+b)/2;elsea=(a+b)/2;enditer=iter+1;endx=(a+b)/2;iter=iter+1;求方程的實(shí)根：test9.mclearall;clc;symsxf=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)運(yùn)行結(jié)果：[0,pi]區(qū)間的根x1=1.8807； 迭代次數(shù)iter1=[pi,2*pi]區(qū)間的根x2=4.6941； 迭代次數(shù)iter2=20[2*pi,3*pi]區(qū)間的根x3=7.8548；迭代次數(shù)iter3=20[3*pi,4*pi]區(qū)間的根x4=10.9955；迭代次數(shù)iter4=20程序：Newton插值：Newtominter.mfunctionf=Newtominter(x,y,x0)%牛頓插值x插值節(jié)點(diǎn)%y為對(duì)應(yīng)的函數(shù)值%Newtonx_0fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsey的維數(shù)不相等！');return;endf=y1=0;l=1;for(i=1:n-1)for(j=i+1:n)y1(j)=(y(j)-y(i))/(x(j)-x(i));endc(i)=y1(i+1);l=l*(t-x(i));f=f+c(i)*l;simplify(f);y=if(i==n-1)if(nargin3) 3個(gè)參數(shù)則給出插值點(diǎn)的插值結(jié)果f='else%2fcollect(f); 將插值多項(xiàng)式展開(kāi)f=vpa(f,6);endendend用等距節(jié)點(diǎn)做f(x)的Newton插值：test10.mn1=5;n2=10;n3=15;x0=[0:0.01:1];y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距節(jié)點(diǎn),5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace(0,1,n2);%等距節(jié)點(diǎn),10y2=sin(pi.*x2);f02=Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距節(jié)點(diǎn),15y3=sin(pi.*x3);f03=Newtominter(x3,y3,x0);plot(x0,y0,'-r')%原圖holdonplot(x0,f01,'-g')%5個(gè)節(jié)點(diǎn)plot(x0,f02,'-k')%10個(gè)節(jié)點(diǎn)plot(x0,f03,'-b')%15個(gè)節(jié)點(diǎn)legend('原圖','5Newton插值多項(xiàng)式','10Newton插值多項(xiàng)式','15個(gè)Newton插值多項(xiàng)式')運(yùn)行結(jié)果：取不同的節(jié)點(diǎn)做牛頓插值。得到結(jié)果圖像如下：可以看出原圖與插值多項(xiàng)式的圖像近似重合，說(shuō)明插值效果較好。程序：Lagrange插值：Lagrange.mfunction[f,f0]=Lagrange(x,y,x0)%Lagrange插值x為插值結(jié)點(diǎn)，y為對(duì)應(yīng)的函數(shù)值，x0為要計(jì)算的點(diǎn)。%L_n(x)fL_n(x0)f0。symst;if(length(x)==length(y))n=length(x);elsey的維數(shù)不相等！');return;end%檢錯(cuò)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)ll*(t-x(j))/(x(i)-x(j)); Lagrange基函數(shù)end;f=fl; %Lagrange插值函數(shù)simplify(f); %化簡(jiǎn)if(i==n)if(nargin==3)0='; 如果3個(gè)參數(shù)則計(jì)算插值點(diǎn)的函數(shù)值elsefcollect(f); %2個(gè)參數(shù)則將插值多項(xiàng)式展開(kāi)fvpa(f,6); %6位精度的小數(shù)endendend用等距節(jié)點(diǎn)做Lagrange插值：test11.mclearall;clc;n1=5;n2=10;n3=15;x0=[-5:0.02:5];y0=1./(1+x0.^2);x1=linspace(-5,5,n1);%等距節(jié)點(diǎn),5y1=1./(1+x1.^2);[f1,f01]=Lagrange(x1,y1,x0);x2=linspace(-5,5,n2);%等距節(jié)點(diǎn),10y2=1./(1+x2.^2);[f2,f02]=Lagrange(x2,y2,x0);x3=linspace(-5,5,n3);%等距節(jié)點(diǎn),15y3=1./(1+x3.^2);[f3,f03]=Lagrange(x3,y3,x0);plot(x0,y0,'-r')%原圖holdonplot(x0,f01,'-b')%5個(gè)節(jié)點(diǎn)plot(x0,f02,'-g')%10個(gè)節(jié)點(diǎn)plot(x0,f03,'-k')%15個(gè)節(jié)點(diǎn)xlabel('x');ylabel('f(x),L(x)');legend('原圖','5Lagrange插值多項(xiàng)式','10Lagrange插值多項(xiàng)式','15Lagrange插值多項(xiàng)式')運(yùn)行結(jié)果：選取了5,10,15個(gè)節(jié)點(diǎn)做Lagrange插值，得到原圖與插值多項(xiàng)式的圖像如下：當(dāng)選取節(jié)點(diǎn)數(shù)較多時(shí)，選取15個(gè)節(jié)點(diǎn)時(shí)出現(xiàn)Runge現(xiàn)象。程序：復(fù)合梯形求積：Comtrap.mfunctions=Comtrap(f,a,b,n)%復(fù)合梯形求積s積分結(jié)果%f積分函數(shù)%a，b積分區(qū)間%n區(qū)間個(gè)數(shù)h=(b-a)/n;index=(a+h):h:(b-h);s1=sum(subs(f,index));s=(h/2)*(subs(f,a)+2*s1+subs(f,b));復(fù)合Simpson求積：functions=Simpson(f,a,b,n)%Simpson求積s積分結(jié)果%f積分函數(shù)%a，b積分區(qū)間%n區(qū)間個(gè)數(shù)h=(b-a)/(2*n);index1=(a+h):(2*h):(b-h);index2=(a+2*h):(2*h):(b-2*h);s1=sum(subs(f,index1));s2=sum(subs(f,index2));s=(h/3)*(subs(f,a)+4*s1+2*s2+subs(f,b));計(jì)算f(x)積分：test12.mclearall;clc;symsxf=exp(3.*x).*cos(pi.*x);a=0;b=2*pi;disp('精確積分值');I=vpa(int(f,x,a,b),10)n1=50;n2=100;n3=200;n4=500;n5=1000;disp('50,100,200,500,1000的復(fù)合梯形積分值');T1=vpa(Comtrap(f,a,b,n1),10)et1=abs(I-T1)T2=vpa(Comtrap(f,a,b,n2),10)et2=abs(I-T2)T3=vpa(Comtrap(f,a,b,n3),10)et3=abs(I-T3)T4=vpa(Comtrap(f,a,b,n4),10)et4=abs(I-T4)T5=vpa(Comtrap(f,a,b,n5),10)et5=abs(I-T5)disp('50,100,200,500,1000Simpson積分值');s1=vpa(Simpson(f,a,b,n1),10)es1=abs(I-s1)s2=vpa(Simpson(f,a,b,n2),10)es2=abs(I-s2)s3=vpa(Simpson(f,a,b,n3),10)es3=abs(I-s3)s4=vpa(Simpson(f,a,b,n4),10)es4=abs(I-s4)s5=vpa(Simpson(f,a,b,n5),10)es5=abs(I-s5)運(yùn)行結(jié)果：精確積分值I=35232483.36復(fù)合梯形與復(fù)合Simpson求積與誤差區(qū)間個(gè)數(shù)n復(fù)合梯形求積T誤差eT復(fù)合Simpson求積誤差eS5035125341.19107142.1735231407.871075.4910035204891.227592.1635232416.2467.1220035225534.986948.3835232479.174.1950035231369.371113.9935232483.250.11100035232204.78278.5835232483.360.0程序：GaussLegendre求積：GaussLegendre.mfunctions=GaussLegendre(f,a,b,n)%--Gauss-Legendre2-5個(gè)求積結(jié)點(diǎn)%f 被積函數(shù)%b,a積分上下限%n n+1%s 積分結(jié)果ta=(b-a)/2;tb=(a+b)/2;switchncase1,subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb));case2,s=ta*(0.5555556*subs(sym(f),findsym(sym(f)),ta*0.7745967+tb)+...0.8888889*subs(sym(f),findsym(sym(f)),tb));case3,s=ta*(0.3478548*subs(sym(f),findsym(sym(f)),ta*0.8611363+tb)+...0.6521452*subs(sym(f),findsym(sym(f)),-ta*0.3399810+tb));case4,s=ta*(0.2369269*subs(sym(f),findsym(sym(f)),ta*0.9061798+tb)+...0.5688889*subs(sym(f),findsym(sym(f)),tb));endGaussChebyshev求積：GaussChebyshev.mfunctions=GaussChebyshev(f,n)%GaussChebyshev求積s積分值%f積分函數(shù)%求積結(jié)點(diǎn)個(gè)數(shù)n+1k=0:n;x=cos((2.*k+1).*pi/(2.*(n+1)));a=pi/(n+1);s=a*sum(subs(f,x));f(x)積分：test13.mclearall;clc;symsxf1=x.^2;f2=sin(x)./x;a=0;b=pi/2;disp('第一個(gè)實(shí)際積分值');I1=vpa(int(f1/p,x,-1,1),8)disp('第二個(gè)實(shí)際積分值');I2=vpa(int(f2,x,a,b),8)disp('2,3,5Guass型積分');disp('2,3,5Guass型積分');GL2=vpa(GaussLegendre(f2,a,b,n1-1),8)GL3=vpa(GaussLegendre(f2,a,b,n2-1),8)GL5=vpa(GaussLegendre(f2,a,b,n3-1),8)運(yùn)行結(jié)果：（1）實(shí)際積分值：I1=1.57079632,3,5點(diǎn)GaussChebyshev型積分I2=1.5707963I3=1.5707963I5=1.5707963（2）實(shí)際積分值I2=1.37076222,3,5點(diǎn)GaussLegendreGL2=1.370419GL3=1.3707635GL5=1.3707622程序：Euler法：Euler.mfunctionT=Euler(f,x0,xn,y0,h)%Euler法%x,y微分方程的解%f微分方程%x0,y0初始值%xn端點(diǎn)%h步長(zhǎng)n=(xn-x0)/h;%區(qū)間個(gè)數(shù)x=zeros(1,n+1);x(1)=x0;y(1)=y0;fork=1:ny(k+1)=y(k)+h*feval(f,x(k),y(k));x(k+1)=x(k)+h;endT=[x',y'];改進(jìn)Euler法：TranEuler.m