《数值分析报告编程及运行结果.docx》由会员分享，可在线阅读，更多相关《数值分析报告编程及运行结果.docx（30页珍藏版）》请在三一办公上搜索。
1、高斯消元法1.程序：clearformat ratA=input(输入增广矩阵A=)m,n=size(A);for i=1:(m-1)numb=int2str(i);disp(第,numb,次消元后的增广矩阵)for j=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endAend%回代过程disp(回代求解)x(m)=A(m,n)/A(m,m);for i=(m-1):-1:1x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)/A(i,i);end2.运行结果:高斯选列主元消元法1.程序：clearformat ratA=input(输入增
2、广矩阵A=)m,n=size(A);for i=1:(m-1)numb=int2str(i);disp(第,numb,次选列主元后的增广矩阵)temp=max(abs(A(i:m,i);a,b=find(abs(A(i:m,i)=temp);tempo=A(a(1)+i-1,:);A(a(1)+i-1,:)=A(i,:);A(i,:)=tempodisp(第,numb,次消元后的增广矩阵)for j=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endAend%回代过程disp(回代求解)x(m)=A(m,n)/A(m,m);x(i)=(A(i,n)-A(
3、i,i+1:m)*x(i+1:m)/A(i,i);endx2.运行结果:# Command WindowXFile Edit DebugDeslctop WindowHelpB .输入增广矩降LL 1 6A :-13 1 42-6 1 -5A =1116-13142-61-5第L次迭列主元后的增广拒阵A =2-61 -13141116第L次消元后的增广拒阵1A =2-61 -003/23/2041/217/2第2次迭列主元后的增广拒阵A =2-61 -V精彩文档?OVR追赶法1.程序：function x,L,U=zhuiganfa(a,b,c,f)a=input(输入矩阵-1对角元素a=);
4、b=input(输入矩阵对角元素b=);c=input(输入矩阵+1对角元素c=);f=input(输入增广矩阵最后一列元素f=);n=length(b);%对A进行分解u(1)=b (1);for i=2:nif(u(i-1)=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);elsebreak;endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;%求解Ly=by(i)=f(i);for i=2:ny(i)=f(i)-l(i-1)*y(i-1);end%求解Ux=yif(u(
5、n)=0)x(n)=y(n)/u(n);end for i=n-1:-1:1x(i)=(y(i)-c (i) *x(i+1)/u(i);end2.运行结果:Command WindowFile Edit Debug DesktopWindowHelp输入拒用-1对角元素莎田2.1 输入拒防对角元素b= 3 X X 3 输入拒阵+L对角元素TL 1.1 输入增广矩阵最后一列元素玲2,1,2, 4 ans15/19-7/1910/1922/19高斯-塞德尔迭代格式1.程序：function x=Gauss_Seidel(a,b)a=input(输入系数矩阵a=)b=input(输入增广矩阵最后一列
6、b=);e=0.5e-7;n=length(b);N=50;x=zeros(n,1);t=zeros(n,1);for k=1:Nsum=0;E=0;t(1:n)=x(1:n);for i=1:nx(i)=(b (i)-a(i,1:(i-1)*x(1:(i-1)-a(i,(i+1):n)*t(i+1):n)/a(i,i);endif norm(x-t)ekbreak;endend2.运行结果:雅戈比迭代格式1.程序：function x=Jocabi(a,b)a=input(输入系数矩阵a=);b=input(输入增广矩阵最后一列b=);e=0.5e-7;n=length(b);N=100;x
7、=zeros(n,1);y=zeros(n,1);for k=1:Nsum=0;for i=1:ny(i)=(b(i)-a(i,1:n)*x(1:n)+a(i,i)*x (i) )/a(i,i);endfor i=1:nsum=sum+(y(i)-x(i)A2;endif sqrt(sum)ekbreak;elsefor i=1:nx(i)=y (i)；endendendif k=N warning(未能找到近似解);end2.运行结果：逐次超松弛法（SOR）1.程序：function n,x=sor22（A,b,X,nm,w,ww）%用超松弛迭代法求解方程组Ax=b%输入：A为方程组的系数矩
8、阵，b为方程组右端的列向量，X为迭 代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为 松弛因子%输出：x为求得的方程组的解构成的列向量，n为迭代次数 A=input（输入系数矩阵A=）;b=input（输入方程组右端的列向量b=）;X=input（输入迭代初值构成的列向量X=）;nm=input（输入最大迭代次数nm=）;w=input（输入误差精度w=）;ww=input（输入松弛因子ww=）;n=1;m=length（A）;D=diag(diag(A);%令A=D-L-U,计算矩阵DL=tril(-A)+D;%令A=D-L-U,计算矩阵LU=triu(-A)+D;%令A=D-L
9、-U,计算矩阵UM=inv(D-ww*L)*(1-ww)*D+ww*U);% 计算迭代矩阵g=ww*inv(D-ww*L)*b;%计算迭代格式中的常数项%下面是迭代过程while n=nmx=M*X+g;%用迭代格式进行迭代if norm（x-X,inf）epsx=(a+b)/2;%检查是否大于值if (xA3)-3*x-1)Mb=xelsea=xendk=k+1end2.运行结果:Newton迭代法(切线法)1.程序：function x=nanewton(fname,dfname,x0,e,N)%newton迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内
10、嵌函数，x0为迭代初值，e为精度要求x=x0;x0=x+2*e;k=0;if nargin5,N=500;endif nargine&kA VLAcA * A*LUJ *4 X * V L4 XU A * !* 4 A4ULMJ * 4 X*,_A fun=iniine( x-cos(x) );dfun=iniine( 1+sin(x); K=nanevton(furL? dfun, 1, le-5, 100)E ? Undef ined funct ion or method , nanewton, for input arguments of type K=naneTton(fun; df
11、un, 1, le-5, 10j)0.75040.73910.73910.73910.7391v|OVR .：割线方式迭代法1.程序：function x=ge_xian_fa(fname,dfname,x0,x1,e,N)%割线方式迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数，x0,x1分别为迭代初值，e为精度要求k=0;a=x1;b=x0;if nargin5,N=500;endif nargine&k0,x0=x;b=x0;elsex1=x;a=x1;endx=x1-(x1-x0)/(feval(fname,x1)-feval(fname,x0
12、)*feval(fname,x1);numb=int2str(k);disp(第,numb,次计算后 x=)fprintf(%fnn,x);endif k=N,warning（已达迭代次数上限）;end2.运行结果:Newton插值1.程序：%保存文件名为New_Int.m%Newton基本插值公式%x为向量，全部的插值节点%y为向量，差值节点处的函数值%xi为标量，是自变量%yi为xi出的函数估计值function yi=newton_chazhi(x,y,xi)n=length(x);m=length(y);if n=merror(The lengths of X ang Y must b
13、e equal!);return;end%计算均差表YY=zeros(n);Y(：,1)=y;for k=1:n-1for i=1:n-kif abs(x(i+k)-x(i) clear yi=nevton_chazhi(Lpi/6:pi/41? ZO. 5;0. 7071, 2*pi/9)J ? Undef ined funct ion or method ? newton_chazhiJ for input arguments of ty _ yi=newton_chazhi(pi/6:pi/4Z, Z0. 5;0. 7071, 2*pi/9)yi =0.6381 yi=nevton_ch
14、azhi(pi/6:pi/4;pi/31, 10.5;0. 7071;0.86601,2*pi/9)yi =0. 6434Lagrange寸插值1.程序：function y0 = Language(x,y,x0)syms t l;if length(x)=length(y)n = length(x);elsedisp(x和y的维数不相等！)；return;% 检错endh=sym(0);for i=1:nl=sym(y (i);for j=1:i-1l=l*(t-x(j)/(x(i)-x (j);end;for j=i+1:nl=l*(t-x(j)/(x(i)-x (j);end;h=h+l
15、;endsimplify(h);if nargin = 3y0 = subs (h,t,x0); %计算插值点的函数值elsey0=collect(h);y0 = vpa(y0,6); %将插值多项式的系数化成6位精度的小数 end2.运行结果:或 Command WindowCZ1回*File Edit Debug Desktop Window 旦elp叔-七-11,1, 6； yO=Language(% y)yO 二 500000*t3+2. *t+2. SOOOO+t+l.最小二乘法1.程序：function p=nafit(x,y,m)%多项式拟合% x, y为已知数据点向量，分别表示横，纵坐标，m为拟合多项式的 次数，结果返回m次拟合多项式系数，从高次到低次存放在向量p 中.A=zeros(m+1,m+1);for i=0:mfor j=0:mA(i+1,j+1)=sum(x.A(i+j)；endb(i+1)=sum(x.Ai.*y)；enda=Ab;p=fliplr(a);t=0:0.1:1.6;S3=polyval(p,t);plot(x,y,p k);hold onplot(t,S3,r);2.运行结果:
