MATLAB实验报告整合.docx
- 文档编号:23538869
- 上传时间:2023-05-18
- 格式:DOCX
- 页数:21
- 大小:172.19KB
MATLAB实验报告整合.docx
《MATLAB实验报告整合.docx》由会员分享,可在线阅读,更多相关《MATLAB实验报告整合.docx(21页珍藏版)》请在冰豆网上搜索。
MATLAB实验报告整合
数值分析
MATLAB实验报告
姓名:
张宗豪
学号:
1304100514
指导老师:
易昆南
专业班级:
统计1003
具体内容:
学号
1304100514
班级
统计1003
姓名
张宗豪
指导教师
易昆南
实验题目1
用Gauss消元法、Gauss列主元消去法、LU分解法(Doolittle法)解方程组
评分
1、设计(实习)目的:
1.了解数值分析在实际问题中的应用
2.通过实践加深对这门课程的了解
3.熟悉简单程序结构,如循环结构(for循环、while循环)选择结构(if-else-if)、分支语句(switch-case-otherwise)。
2、实验内容:
求解方程组:
2.5x1+2.3x2-5.1x3=3.7
5.3x1+9.6x2+1.5x3=3.8
8.1x1+1.7x2-4.3x3=5.5
3.详细设计:
(一)Gauss消元法:
functionx=Gauss(A,b)
n=length(b);
fori=1:
n-1
fork=i+1:
n
forj=i+1:
n
ifabs(A(i,i))<10^(-6)
warning('分母为零不能计算!
');
return;
elseA(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);
end
end
b(k)=b(k)-b(i)*A(k,i)/A(i,i);
A(k,i)=0;
end
end
x(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A(i,j)*x(j);
end
x(i)=(b(i)-sum)/A(i,i);
end
(二)Gauss列主元消去法:
functionx=Gaussrow(A,b)
n=length(b);
x=zeros(n,1);
c=zeros(1,n);
t=0;
fori=1:
n-1
max=abs(A(i,i));
m=i;
forj=i+1:
n
ifmax max=abs(A(j,i)); m=j; end end ifm~=i fork=1: n c(k)=A(i,k); A(i,k)=A(m,k); A(m,k)=c(k); end t=b(i); b(i)=b(m); b(m)=t; end fork=i+1: n forj=i+1: n A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); end b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0; end end x(n)=b(n)/A(n,n); fori=n-1: -1: 1 sum=0; forj=i+1: n sum=sum+A(i,j)*x(j); end x(i)=(b(i)-sum)/A(i,i); end (三)LU分解法(Doolittle法): function[x,L,U]=Doolittle(A,b) %LU分解 n=length(b); U=zeros(n,n); L=eye(n,n); U(1,: )=A(1,: ); L(2: n,1)=A(2: n,1)/U(1,1); fork=2: n U(k,k: n)=A(k,k: n)-L(k,1: k-1)*U(1: k-1,k: n); L(k+1: n,k)=(A(k+1: n,k)-L(k+1: n,1: k-1)*U(1: k-1,k))/U(k,k); end %解下三角方程组Ly=b y=zeros(n,1); y (1)=b (1); fork=2: n y(k)=b(k)-L(k,1: k-1)*y(1: k-1); end x=zeros(n,1); x(n)=y(n)/U(n,n); fork=n-1: -1: 1 x(k)=(y(k)-U(k,k+1: n)*x(k+1: n))/U(k,k); end 4: 实验结果: 输入A=[2.5,2.3,-5.1;5.3,9.6,1.5;8.1,1.7,-4.3]; b=[3.7,3.8,5.5]'; (一)Gauss消元法: >>x=Gauss(A,b) x= 0.40670.2368-0.4193 (二)Gauss列主元消去法: >>x=Gaussrow(A,b) x= 0.4067 0.2368 -0.4193 (三)LU分解法(Doolittle法): >>[x,L,U]=Doolittle(A,b) x= 0.4067 0.2368 -0.4193 L= 1.000000 2.12001.00000 3.2400-1.21761.0000 U= 2.50002.3000-5.1000 04.724012.3120 0027.2152 学号 1304100514 班级 统计1003 姓名 张宗豪 指导教师 易昆南 实验题目2 用牛顿法,复化梯形法,复化辛普森法,复化辛普森变步长法,求积分 评分 1、设计(实习)目的: 1.了解数值分析在实际问题中的应用 2.通过实践加深对这门课程的了解 3.熟悉简单程序结构,如循环结构(for循环、while循环)选择结构(if-else-if)、分支语句(switch-case-otherwise)。 2、实验内容: 用牛顿法求6*x.^3+5*x.^2+4*x+3的积分并画图,用复化梯形法,复化辛普森法,复化辛普森变步长法求lg(1+x)的积分 3.详细设计: 一.牛顿法: functiony=newton(a,b,n) x=a: (b-a)/n: b;%插值节 y=6*x.^3+5*x.^2+4*x+3; plot(x,y,'b')%用蓝色线作被插函数图象 holdon n=length(x); z=a: (b-a)/(2*n): b; forj=2: n fori=n: -1: j y(i)=(y(i)-y(i-1))/(x(i)-x(i-j+1));%计算差商 end end u=y(n); m=length(z); forj=1: m fori=n-1: -1: 1 u=y(i)+u*(z(j)-x(i));%计算牛顿插值多项式的值 v(j)=u; end u=y(n); end plot(z,v,'r')%用红色线作牛顿插值多项式图象 holdon 二.复化梯形法: functionTn=comtr(a,b,n) y=inline('log(1+x)','x'); h=(b-a)/n; fork=0: n x(k+1)=a+k*h; ifx(k+1)==0 x(k+1)=10^(-10); end end T1=h/2*(y(x (1))+y(x(n+1))); fori=2: n F(i)=h*y(x(i)); end T2=sum(F); Tn=T1+T2; 三.复化辛普森法: functionSn=comsim(a,b,n) y=inline('log(1+x)','x'); h=(b-a)/n; fork=0: n; x(k+1)=a+k*h; x_k(k+1)=x(k+1)+1/2*h; if(x(k+1)==0)|(x_k(k+1)==0) x(k+1)=10^(-10); x_k(k+1)=10^(-10); end end S1=h/6*(y(x (1))+y(x(n+1))); fori=2: n F1(i)=h/3*y(x(i)); end forj=1: n F2(j)=2*h/3*y(x_k(j)); end S2=sum(F1)+sum(F2); Sn=S1+S2; 四.复化辛普森变步长法: function[Q,cnt]=quad(funfcn,a,b,tol,trace)%funfcn用内置函数输入 ifnargin<4,tol=1.e-3;trace=0;end ifnargin<5,trace=0;end c=(a+b)/2; x=[abca: (b-a)/10: b]; y=feval(funfcn,x); fa=y (1); fb=y (2); fc=y(3); iftrace lims=[min(x)max(x)min(y)max(y)]; ifany(imag(y)) lims(3)=min(min(real(y)),min(imag(y))); lims(4)=max(max(real(y)),max(imag(y))); end ind=find(isfinite(lims)); ifisempty(ind) [mind,nimd]=size(ind); Lims(ind)=1.e30*(-ones(mind,nind).^rem(ind,2)); end axis(lims); plot([c,b],[fc,fb],'.') holdon ifany(imag(y)) plot([c,b],imag([fc,fb]),'+') end end lev=1; ifany(imag(y)) Q0=le30; else Q0=inf; end [Q,cnt]=quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace); cnt=cnt+3; iftrace holdoff axis('auto'); end function[Q,cnt]=quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace) LEVMAX=10; iflev>LEVMAX disp('Recursionlevellimitreachedinquad.Singularitylikely.') Q=Q0; Cnt=0; c=(a+b)/2; iftrace yc=feval(funfcn,c); plot(c,yc,'*'); ifany(imag(yc)) plot(c,imag(yc),'*'); end end else h=b-a; c=(a+b)/2; x=[a+h/4,b-h/4]; f=feval(funfcn,x); iftrace plot(x,f,'.'); ifany(imag(f)) plot(x,imag(f),'+') end end cnt=2; Q1=h*(fa+4*f (1)+fc)/12; Q2=h*(fc+4*f (2)+fb)/12; Q=Q1+Q2; ifabs(Q-Q0)>tol*abs(Q) [Q1,cnt1]=quadstp(funfcn,a,c,tol/2,lev+1,fa,f (1),fc,Q1,trace); [Q2,cnt2]=quadstp(funfcn,c,b,tol/2,lev+1,fc,f (2),fb,Q2,trace); Q=Q1+Q2; Cnt=cnt+cnt1+cnt2; end end 4: 实验结果: 一.牛顿法: 输入newton(0,1,10),得到 ans= Columns1through8 3.00004.56006.80006.00000.0000-0.00000.0000-0.0000 Columns9through11 0.00000.0000-0.0000 二.复化梯形法: 输入comtr(0,1,10),得到 ans= 0.3859 三.复化辛普森法: 输入comsim(0,1,10),得到 ans= 0.3830 四.复化辛普森变步长法: 输入quad(inline('log(1+x)','x'),0,1,10),得到 ans= 0.3863 学号 1304100514 班级 统计1003 姓名 张宗豪 指导教师 易昆南 实验题目3 求一阶微分f=(x-y)/2x=[03]y(0)=1区间间隔分别为h=11/21/41/8的值。 评分 1、设计(实习)目的: 1.了解数值分析在实际问题中的应用 2.通过实践加深对这门课程的了解 3.熟悉简单程序结构,如循环结构(for循环、while循环)选择结构(if-else-if)、分支语句(switch-case-otherwise)。 2、实验内容: 求一阶微分f=(x-y)/2,x=[03],y(0)=1,区间间隔分别为h=1,1/2,1/4,1/8的值。 3.详细设计: 一阶微分m.文件: functionf=myfun(x,y) f=(x-y)/2; 一.欧拉法: function[xout,yout]=euler(ypfun,xspan,y0,h) f=ypfun; x=xspan (1): h: xspan (2); n=length(x); y=zeros(1,n); y (1)=y0; fori=1: n-1 y(i+1)=y(i)+h*feval(f,x(i),y(i)); end xout=x; yout=y; plot(x,y,'-r') 二.修恩法(改进的欧拉法): function[xout,yout]=eulerpro(ypfun,xspan,y0,h) f=ypfun; x=xspan (1): h: xspan (2); n=length(x); y=zeros(1,n); y (1)=y0; fori=1: n-1 yp=y(i)+h*feval(f,x(i),y(i)); yc=y(i)+h*feval(f,x(i+1),yp); y(i+1)=(yp+yc)/2; end xout=x; yout=y; plot(x,y,'-r') 三.四阶泰勒展开式法: function[xout,yout]=taylor4(xspan,y0,h) f=inline('(x-y)/2'); f1=inline('(2-x+y)/4'); f2=inline('(-2+x-y)/8'); f3=inline('(2-x+y)/16'); x=xspan (1): h: xspan (2); n=length(x); y=zeros(1,n); y (1)=y0; fori=1: n-1 y(i+1)=y(i)+h*f(x(i),y(i))+h.^2/2*f1(x(i),y(i))+h.^3/6*f2(x(i),y(i))+... h.^4/24*f3(x(i),y(i)); end xout=x; yout=y; plot(x,y,'-r') 四.四阶龙格-库特法: function[xout,yout]=varargout4(ypfun,xspan,y0,h) f=ypfun; x=xspan (1): h: xspan (2); n=length(x); y=zeros(1,n); y (1)=y0; fori=1: n-1 f1=f(x(i),y(i)); f2=f(x(i)+h/2,y(i)+h/2*f1); f3=f(x(i)+h/2,y(i)+h/2*f2); f4=f(x(i)+h,y(i)+h*f3); y(i+1)=y(i)+h/6*(f1+2*f2+2*f3+f4); end xout=x; yout=y; plot(x,y,'-r') 五.精确解: function[xout,yout]=exact(xspan,y0,h) f=inline('3*exp(-x/2)-2+x'); x=xspan (1): h: xspan (2); n=length(x); y (1)=y0; fori=1: n-1 y(i+1)=f(x(i+1)); end xout=x; yout=y; plot(x,y,'-r') 4: 实验结果: 一.欧拉法: 输入[x,y]=euler('myfun',[03],1,1) x= 0123 y= 1.00000.50000.75001.3750 二.修恩法(改进的欧拉法): 输入[x,y]=eulerpro('myfun',[03],1,1) x= 0123 y= 1.00000.87501.17191.7324 三.四阶泰勒展开式法: 输入[x,y]=taylor4([03],1,1) x= 0123 y= 1.00000.82031.10451.6702 四.四阶龙格-库特法: 输入[x,y]=varargout4('myfun',[03],1,1) x= 0123 y= 1.00000.82031.10451.6702 五.精确解: 输入[x,y]=exact([03],1,1) x= 0123 y= 1.00000.81961.10361.6694 当h=1/2,1/4,1/8时,同理可得,此处略。 5: 实验总结 通过实验了解到,数值分析是一门理论和实际并重的基础学科,掌握好数值分析的基本知识,对以后的学习和工作都有很大的帮助和指导意义。 我觉得收获呢,收获了一些基础性的数值分析知识,对于以后的学习很有帮助,感觉很高兴。 特别是解线性方程组,实用价值很高;另外,在实验的过程中,我发现复化辛普森法和复化辛普森变步长法要比其他方法精度更高一些;了解了除Matlab自身求微分函数以外的数值分析中的其他求微分的方法,实用性很强,这也是一种收获吧。 我立志要学好这门课程。 虽然现在我还不完全熟悉这门课程,但我对自己有信心,相信我能做到。 姓名: 张宗豪 2012年04月21日
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 实验 报告 整合