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

矩阵与数值分析上机作业学院：班级：姓名：学号：授课老师：注：编程语言Matlab程序：Norm.m函数functions=Norm(x,m)%x的范数%m1,2,inf1,2，无穷范数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%无穷-范数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时');disp('x1-范数:');disp(Norm(x1,1));disp('x2-范数:');disp(Norm(x1,2));disp('x的无穷-范数:');disp(Norm(x1,inf));disp('y1-范数:');disp(Norm(y1,1));disp('y2-范数:');disp(Norm(y1,2));disp('y的无穷-范数:');disp(Norm(y1,inf));disp('n=100时');disp('x1-范数:');disp(Norm(x2,1));disp('x2-范数:');disp(Norm(x2,2));disp('x的无穷-范数:');disp(Norm(x2,inf));disp('y1-范数:');disp(Norm(y2,1));disp('y2-范数:');disp(Norm(y2,2));disp('y的无穷-范数:');disp(Norm(y2,inf));disp('n=1000时');disp('x1-范数:');disp(Norm(x3,1));disp('x2-范数:');disp(Norm(x3,2));disp('x的无穷-范数:');disp(Norm(x3,inf));disp('y1-范数:');disp(Norm(y3,1));disp('y2-范数:');disp(Norm(y3,2));disp('y的无穷-范数:');disp(Norm(y3,inf));运行结果：n=10时x的1-范:2.9290；x的2-范:1.2449；x的无-范:1y的1-范:55； y的2-范:19.6214；y的无-范n=100时x1:5.1874；x2-范数:1.2787；x:1y1:5050；y2-:581.6786；y:100n=1000时x的1-范:7.4855；x的2-范数:1.2822； x的无范:1y1500500；y2-范数:1.8271e+004；y:1000程序Test2.mclearclc;h=2*10^(-15)/n;%步长x=-10^(-15):h:10^(-15);%第一种原函数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('第二种算法');运行结果：显然第二种算法结果不准确计算机进行舍入造成d11。程序：秦九韶算法:QinJS.mfunctiony=QinJS(a,x)%y输出函数值%a多项式系数，由高次到零次%x给定点n=length(a);s=a(1);fori=2:ns=s*x+a(i);endy=s;计算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)');运行结果：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)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('当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的逆矩阵的准确逆矩阵:');InvA=inv(A)InvAL=zeros(n)LUA的逆矩阵I=eye(n);fori=1:nInvAL(:,i)=U\(L\I(:,i));endLUA的逆矩阵:');InvALEnd运行结果：n=5,6,7当n=5方程精确解:x=1.00000.50000.33330.25000.2000LUxLU=1.00000.50000.33330.25000.2000PLUxPLU=1.00000.50000.33330.25000.2000当n=6方程精确解:x=1.00000.50000.33330.25000.20000.1667LUxLU=1.00000.50000.33330.25000.20000.1667PLUxPLU=1.00000.50000.33330.25000.20000.1667当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当n=5A的准确逆矩阵: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当n=6A的准确逆矩阵: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当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.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计算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运行结果：只列出了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)运行结果：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%--迭代步数控制M%--·返回求得的解x%--·返回迭代步数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=D\b;x=J*x0+fn=1;%迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=J*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M):迭代次数太多，可能不收敛?'return;endendGauss_Seidel迭代：Gauss_Seidel.mfunction[x,n]=Gauss_Seidel(A,b,x0)%--Gauss-Seidel迭代法解线性方程组%--方程组系数阵A%--方程组右端项b%--初始值x0%--求解要求的精确度eps%--迭代步数控制M%--返回求得的解x%--返回迭代步数neps=1.0e-5;M=10000;D=diag(diag(A));A的对角矩阵L=-tril(A,-1);A的下三角阵U=-triu(A,1);A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+fn=1; 迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=G*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M):迭代次数太多，可能不收敛！'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迭代过程：');[xj,nj]=Jaccobi(A,b,x0);disp('Jacobi最终迭代结果：');xjdisp('迭代次数');njdisp('Gauss-Seidel迭代过程：');[xg,ng]=Gauss_Seidel(A,b,x0);disp('Gauss-Seidel最终迭代结果：');xgdisp('迭代次数');ng运行结果：精确解x=-0.0820-1.80331.1311迭代初始值x0=000Jacobi迭代过程：x=-0.40000.50000.6667err=0.6667x=0.1000-1.03330.4533err=1.5333...x=-0.0820-1.80331.1311err=9.6603e-006Jacobi最终迭代结果：xj=-0.0820-1.80331.1311迭代次数nj=281Gauss-Seidel迭代过程：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最终迭代结果：xg=-0.0820-1.80331.1311迭代次数ng=20程序：Newton迭代法：Newtoniter.mfunction[x,iter,fvalue]=Newtoniter(f,df,x0,eps,maxiter)%牛顿法x得到的近似解%iter迭代次数%fvaluex处的值%f，df被求的非线性方程及导函数%x0初始值%eps允许误差限%maxiter最大迭代次数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迭代次数%fvaluex处的值%f被求的非线性方程%x0，x1初始值%eps允许误差限%maxiter最大迭代次数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求方程的实根：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迭代结果');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)disp('弦截法初始值');xk2_0=5xk2_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.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弦截法迭代结果：取两种不同初值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迭代次数%f求解的方程%a，b求解区间%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;求方程的实根：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)运行结果：[0,pi]区间的根x1=1.8807； 迭代次数iter1=[pi,2*pi]区间的根x2=4.6941； 迭代次数iter2=20[2*pi,3*pi]区间的根x3=7.8548；迭代次数iter3=20[3*pi,4*pi]区间的根x4=10.9955；迭代次数iter4=20程序：Newton插值：Newtominter.mfunctionf=Newtominter(x,y,x0)%牛顿插值x插值节点%y为对应的函数值%Newtonx_0fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsey的维数不相等！');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个参数则给出插值点的插值结果f='else%2fcollect(f); 将插值多项式展开f=vpa(f,6);endendend用等距节点做f(x)的Newton插值：test10.mn1=5;n2=10;n3=15;x0=[0:0.01:1];y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距节点,5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace(0,1,n2);%等距节点,10y2=sin(pi.*x2);f02=Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距节点,15y3=sin(pi.*x3);f03=Newtominter(x3,y3,x0);plot(x0,y0,'-r')%原图holdonplot(x0,f01,'-g')%5个节点plot(x0,f02,'-k')%10个节点plot(x0,f03,'-b')%15个节点legend('原图','5Newton插值多项式','10Newton插值多项式','15个Newton插值多项式')运行结果：取不同的节点做牛顿插值。得到结果图像如下：可以看出原图与插值多项式的图像近似重合，说明插值效果较好。程序：Lagrange插值：Lagrange.mfunction[f,f0]=Lagrange(x,y,x0)%Lagrange插值x为插值结点，y为对应的函数值，x0为要计算的点。%L_n(x)fL_n(x0)f0。symst;if(length(x)==length(y))n=length(x);elsey的维数不相等！');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)ll*(t-x(j))/(x(i)-x(j)); Lagrange基函数end;f=fl; %Lagrange插值函数simplify(f); %化简if(i==n)if(nargin==3)0='; 如果3个参数则计算插值点的函数值elsefcollect(f); %2个参数则将插值多项式展开fvpa(f,6); %6位精度的小数endendend用等距节点做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);%等距节点,5y1=1./(1+x1.^2);[f1,f01]=Lagrange(x1,y1,x0);x2=linspace(-5,5,n2);%等距节点,10y2=1./(1+x2.^2);[f2,f02]=Lagrange(x2,y2,x0);x3=linspace(-5,5,n3);%等距节点,15y3=1./(1+x3.^2);[f3,f03]=Lagrange(x3,y3,x0);plot(x0,y0,'-r')%原图holdonplot(x0,f01,'-b')%5个节点plot(x0,f02,'-g')%10个节点plot(x0,f03,'-k')%15个节点xlabel('x');ylabel('f(x),L(x)');legend('原图','5Lagrange插值多项式','10Lagrange插值多项式','15Lagrange插值多项式')运行结果：选取了5,10,15个节点做Lagrange插值，得到原图与插值多项式的图像如下：当选取节点数较多时，选取15个节点时出现Runge现象。程序：复合梯形求积：Comtrap.mfunctions=Comtrap(f,a,b,n)%复合梯形求积s积分结果%f积分函数%a，b积分区间%n区间个数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));复合Simpson求积：functions=Simpson(f,a,b,n)%Simpson求积s积分结果%f积分函数%a，b积分区间%n区间个数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));计算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的复合梯形积分值');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)运行结果：精确积分值I=35232483.36复合梯形与复合Simpson求积与误差区间个数n复合梯形求积T误差eT复合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个求积结点%f 被积函数%b,a积分上下限%n n+1%s 积分结果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积分函数%求积结点个数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('第一个实际积分值');I1=vpa(int(f1/p,x,-1,1),8)disp('第二个实际积分值');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)运行结果：（1）实际积分值：I1=1.57079632,3,5点GaussChebyshev型积分I2=1.5707963I3=1.5707963I5=1.5707963（2）实际积分值I2=1.37076222,3,5点GaussLegendreGL2=1.370419GL3=1.3707635GL5=1.3707622程序：Euler法：Euler.mfunctionT=Euler(f,x0,xn,y0,h)%Euler法%x,y微分方程的解%f微分方程%x0,y0初始值%xn端点%h步长n=(xn-x0)/h;%区间个数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'];改进Euler法：TranEuler.m