数值分析课程设计.docx
- 文档编号:22995571
- 上传时间:2023-04-29
- 格式:DOCX
- 页数:17
- 大小:57.90KB
数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享,可在线阅读,更多相关《数值分析课程设计.docx(17页珍藏版)》请在冰豆网上搜索。
数值分析课程设计
1
2用高斯消元法的消元过程作矩阵分解。
设
消元过程可将矩阵A化为上三角矩阵U,试求出消元过程所用的乘数
、
、
并以如下格式构造下三角矩阵L和上三角矩阵U
验证:
矩阵A可以分解为L和U的乘积,即A=LU。
Solution:
function[L,U,flag]=LU_Decom(A)
[n,m]=size(A);
ifn~=m
error('TherowsandcolumnsofmatrxAmustbeequal!
');
return;
end
L=eye(n);U=zeros(n);flag='OK';
fork=1:
n
forj=k:
n
z=0;
forq=1:
k-1
z=z+L(k,q)*U(q,j);
end
U(k,j)=A(k,j)-z;
end
ifabs(U(k,k)) flag='failure';return; end fori=k+1: n z=0; forq=1: k-1 z=z+L(i,q)*U(k,k); end L(i,q)=(A(i,k)-z)/U(k,k); end end 输入命令: >>A=[20,2,3;1,8,1;2,-3,15]; >>[L,U,flag]=LU_Decom(A) 运行结果: L= 1.000000 01.00000 -0.375001.0000 U= 20.00002.00003.0000 08.00001.0000 0016.1250 flag= OK >> 3 4 5 6.用欧拉公式和四阶龙格-库塔法分别求解下列初值问题; (1) 欧拉方法 functionE=Euler_1(fun,x0,y0,xN,N) x=zeros(1,N+1);y=zeros(1,N+1); x (1)=x0;y (1)=y0; h=(xN-x0)/N; forn=1: N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end T=[x',y'] functionz=f(x,y) z=0.9*y/(1+2*x); 输入命令 >>Euler_1('f',0,1,1,10) 运行结果 T= 01.0000 0.10001.0900 0.20001.1718 0.30001.2471 0.40001.3172 0.50001.3831 0.60001.4453 0.70001.5045 0.80001.5609 0.90001.6149 1.00001.6668 >> 四阶龙格-库塔方法 functionR=rk4(f,a,b,ya,N) h=(b-a)/N; T=zeros(1,N+1); Y=zeros(1,N+1); T=a: h: b; Y (1)=ya; forj=1: N k1=h*feval(f,T(j),Y(j)); k2=h*feval(f,T(j)+h/2,Y(j)+k1/2); k3=h*feval(f,T(j)+h/2,Y(j)+k2/2); k4=h*feval(f,T(j)+h,Y(j)+k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6; end R=[T'Y'] functionz=(x,y) z=0.9*y/(1+2*x) 输入命令 >>rk4('f',0,1,1,10) 运行结果 R= 01.0000 0.10001.0855 0.20001.1635 0.30001.2355 0.40001.3028 0.50001.3660 0.60001.4259 0.70001.4828 0.80001.5372 0.90001.5894 1.00001.6395 ans= 01.0000 0.10001.0855 0.20001.1635 0.30001.2355 0.40001.3028 0.50001.3660 0.60001.4259 0.70001.4828 0.80001.5372 0.90001.5894 1.00001.6395 >> (2)欧拉方法 functionE=Euler_1(fun,x0,y0,xN,N) x=zeros(1,N+1);y=zeros(1,N+1); x (1)=x0;y (1)=y0; h=(xN-x0)/N; forn=1: N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end T=[x',y'] functionz=(x,y) z=-x*y/(1+x^2); 输入命令 >>Euler_1('f',0,1,1,10) 运行结果 T= 01.0000 0.10001.0900 0.20001.1718 0.30001.2471 0.40001.3172 0.50001.3831 0.60001.4453 0.70001.5045 0.80001.5609 0.90001.6149 1.00001.6668 >> 四阶龙格-库塔方法 functionR=rk4(f,a,b,ya,N) h=(b-a)/N; T=zeros(1,N+1); Y=zeros(1,N+1); T=a: h: b; Y (1)=ya; forj=1: N k1=h*feval(f,T(j),Y(j)); k2=h*feval(f,T(j)+h/2,Y(j)+k1/2); k3=h*feval(f,T(j)+h/2,Y(j)+k2/2); k4=h*feval(f,T(j)+h,Y(j)+k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6; end R=[T'Y'] functionz=(x,y) z=-x*y/(1+x^2); 输入命令 >>rk4('f',0,1,1,10) 运行结果 R= 01.0000 0.10001.0855 0.20001.1635 0.30001.2355 0.40001.3028 0.50001.3660 0.60001.4259 0.70001.4828 0.80001.5372 0.90001.5894 1.00001.6395 ans= 01.0000 0.10001.0855 0.20001.1635 0.30001.2355 0.40001.3028 0.50001.3660 0.60001.4259 0.70001.4828 0.80001.5372 0.90001.5894 1.00001.6395 >> 7用二分法求下列方程在指定区间[a,b]上的实根近似值: (要求误差不超过0.01) (1)x-sinx-1=0,[a,b]=[1,3]; (2)xsinx=1,[a,b]=[0,1.5]. (1)functionroot=bisect(f,a,b,eps) if(nargin==3) eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end if(f1*f2>0) disp('Á½¶Ëµãº¯ÊýÖµ³Ë»ý´óÓÚ0! '); return; else root=FindRoots(f,a,b,eps); end functionr=FindRoots(f,a,b,eps) f_1=subs(sym(f),findsym(sym(f)),a); f_2=subs(sym(f),findsym(sym(f)),b); mf=subs(sym(f),findsym(sym(f)),(a+b)/2); if(f_1*mf>0) t=(a+b)/2; r=FindRoots(f,t,b,eps); else if(f_1*mf==0) r=(a+b)/2; else if(abs(b-a)<=eps) r=(b+3*a)/4; else s=(a+b)/2; r=FindRoots(f,a,s,eps); end end end 输入命令 >>root=bisect('x-sin(x)-1',1,3) 运行结果 root= 1.9346 (2)functionroot=bisect(f,a,b,eps) if(nargin==3) eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end if(f1*f2>0) disp('Á½¶Ëµãº¯ÊýÖµ³Ë»ý´óÓÚ0! '); return; else root=FindRoots(f,a,b,eps); end functionr=FindRoots(f,a,b,eps) f_1=subs(sym(f),findsym(sym(f)),a); f_2=subs(sym(f),findsym(sym(f)),b); mf=subs(sym(f),findsym(sym(f)),(a+b)/2); if(f_1*mf>0) t=(a+b)/2; r=FindRoots(f,t,b,eps); else if(f_1*mf==0) r=(a+b)/2; else if(abs(b-a)<=eps) r=(b+3*a)/4; else s=(a+b)/2; r=FindRoots(f,a,s,eps); end end end 输入命令 >>root=bisect('x*sin(x)-1',0,1.5) 运行结果 root= 1.1142 8分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。 (1) , (2) , (1)function[x,n]=jacobi(A,b,x0,eps,varargin) ifnargin==3 eps=1.0e-6; M=200; elseifnargin<3 error return elseifnargin==5 M=varargin{1}; end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f; n=1; whilenorm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²! '); end end 输入系数矩阵和右端项 A=[4100000000;1410000000;0000000000;0000000000;0000000000;0000000000;0000000000;0000000000;0000000141;0000000014], A= 4100000000 1410000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000141 0000000014 >>b=[2100000012]' b= 2 1 0 0 0 0 0 0 1 2 输入命令 >>x0=ones(10,1); >>[x,n]=Jacobi(A,b,x0) 输出结果 x= 0.2500 -0.2500 NaN NaN NaN NaN NaN NaN NaN NaN n= 1 (2)function[x,n]=jacobi(A,b,x0,eps,varargin) ifnargin==3 eps=1.0e-6; M=200; elseifnargin<3 error return elseifnargin==5 M=varargin{1}; end D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f; n=1; whilenorm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²! '); end end 输入系数矩阵和右端项 9, 10,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 课程设计