数值分析报告编程及运行结果

1、高斯消元法1.程序:clear format ratA=input('输入增广矩阵A=')m, n=size(A);for i=1:(m-1)nu mb=i nt2str(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)

2、;end2.运行结果:协匚ommand WindowWindowHelp-口XFile Edit DebugDesktop输入増广矩P5A=： 1116-13142-61 -5A -1116-13142-61-5第1次消元后的増广矩阵A -1116042100-S-I'17第2次消元后的增广矩時A =1116042100033回代求網X =3212高斯选列主元消元法1. 程序 :clear format ratA=input(' 输入增广矩阵 A=') m,n=size(A);for i=1:(m-1) numb=int2str(i);disp(' 第 '

3、;,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);for i=(m-1):-1:

4、1x(i)=(A(i, n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); endx2.运行结果：Command Window- XFile Edit Debug Desktop Vindow旦出输入增广拒陆AI 1116A-131412 -61 -5JA =1 116*13L42 -61一 0第】次迭列主元后的増广矩0车A -2 -61_Q-13141 116第1次消元后的増广拒陆1A =2 -61-50 03/23/2C4L/217/2第2次选別主元后的增广矩阵A -2 -61-5OVRJ* Command WindowXFile Edit Debug Desktop

5、Window 旦elp0XXXVA第1次消元后的増广矩阵3示A =2-61L003/23/2041/217/2第2次选列主元后的増广矩阵A =2-61L041/217/2003/23/2第2次消元后的増广矩阵A =2-61041/217/2003/23/2回代求解X 321»|追赶法1. 程序：f=');function x,L,U=zhuiganfa(a,b,c,f) a=input(' 输入矩阵 -1 对角元素 a='); b=input(' 输入矩阵对角元素 b='); c=input(' 输入矩阵 +1 对角元素 c='

6、); f=input(' 输入增广矩阵最后一列元素 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(1)=f(1);for i=2: ny(i)=f(i)-|(i-1)*y(i-1);end%求解Ux=yif(u( n)=0)x( n)=y( n)/u( n);endfor

7、 i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1)/u(i);end2.运行结果：和匚ommand Window-XFile Edit Debug DesktopWindowHelp事輸入矩角元素2, 2.1输入拒阵对角7t«b=3,3. 3,3输入矩眸+1对角元tc=：LLl输入增广矩阵最后一列元素h 2, 4 ans =15/19-7/1910/1922/19»高斯 - 塞德尔迭代格式1. 程序：function x=Gauss_Seidel(a,b)a=input(' 输入系数矩阵 a=')b=input(' 输入增广矩阵最后一

8、列 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;end end2.运行结果: Command Window 口File Edit Debug Desktop Window Help输入系数拒阵尸4, 3,0:3,4, 1：-L 0, 4143034-1

9、-104输入增广矩陈最后一列30,-24124ans =96/31120/31-162/31OVR雅戈比迭代格式1. 程序：function x=Jocabi(a,b)a=input(' 输入系数矩阵 a=');b=input(' 输入增广矩阵最后一列 b=');e=0.5e-7;n=length(b);N=100;x=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

10、)八2;endif sqrt(sum)<ekbreak;elsefor i=1:nx(i)=y(i)；endendendif k=N warning('未能找到近似解');end2.运行结果:林匚ommand WindowFile Edit Debug Desktop Window 旦亡Ipaiis3617/11634289/1108-19163/3667输入系教矩3.0;3,4,-1-1t 0t 4输入增广矩阵最后一刊30.-2432ans =96/31120/31-162/31逐次超松弛法（ SOR）1. 程序：function n,x=sor22(A,b,X,nm,

11、w,ww)%用超松弛迭代法求解方程组 Ax=b %俞入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭 代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松 弛因子%俞出： x 为求得的方程组的解构成的列向量， n 为迭代次数A=input（' 俞入系数矩阵 A='）;b=input（' 俞入方程组右端的列向量 b='）;X=input（' 俞入迭代初值构成的列向量 X='）;nm=input（' 俞入最大迭代次数 nm='）;w=input（' 俞入误差精度 w='）;ww=input('

12、 俞入松弛因子 ww=');n=1;m=length(A);D=diag(diag(A); % L=tril(-A)+D;%U=triu(-A)+D;%令 A=D-L-U, 计算矩阵 D 令 A=D-L-U, 计算矩阵 L 令 A=D-L-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'）<wdisp（' 迭代次数为 '）;nd

13、isp（' 方程组的解为 '）;xreturn;% 上面：达到精度要求就结束程序， 输出迭代次数和方程组的 解endX=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛， 输出警告语句及迭代的最 终结果（并不是方程组的解）disp（' 在最大迭代次数内不收敛 !'）;disp（' 最大迭代次数后的结果为 '）;2. 运行结果：d' Command WindowFile Edit Debug Desktop Window 旦elp »输入系数矩P$A=：4, 3, 0:3, 4,-1 ;-l, 0, 4 输入方程组右端的

14、列向堅=24： 30:-24 输入迭代初值构成的列向1X=1 ;1：1 输入最大迭代次数nm=100输入误差精度w=0. 5e- 输入松弛因子ww=l25 迭代次数为方程组的解为3. 09683.8710-5.2258v二分法求解方程的根1. 程序：%其中a, b表示查找根存在的范围，M表示要求解函数的值%f(x) 表示要求解根的方程%eps表示所允许的误差大小function y=er_fen_fa(a,b,M)k=0;eps=0.05while b-a>epsx=(a+b)/2;% 检查是否大于值if (x八3)-3*x -1)>Mb=xelsea=xendk=k+1 end2

15、.运行结果:f毅 Command WindowFile Edit Debug Desktop Window Help_口XAEirar in => mr fen fa at 6f总val (fnanE, a) : fb=f eval (fnaae, b):? Input argument fname is undefined,Eircr in. => Ei_f£n_fa 肚 6faf eval (fname, a' :fb=f eval fnajue, b):» fun=inline(J x '3-3*x-l ):» x=nabisec

16、t ifun, L 2, 0. Oo)? Undefmed function or method nabisect for input argujnents of t;<pe» xnabisect (fun, L 2, 0, 05)1.9063Newton 迭代法(切线法)1. 程序：function x=nanewton(fname,dfname,x0,e,N)%n ewto n迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数，x0为迭代初值，e为精度要求x=x0;x0=x+2*e;k=0;if nargin<5,N=500;en

17、dif nargin<4 e=1e-4;endwhile abs(x0-x)>e&k<N,k=k+1;x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);disp(x)endif k=N,warning(' 已达迭代次数上限 ');end2.运行结果:矶 Command WindowFile Edit Debug Desktop Jfindow HelpA, P"«S h>> fun-ml me ( xcos 'x.1' ) : dfunFinline (' 1+

18、sm' x ' ):>> x=nanewn oil (fun, dfun, 1, le'5, 100)? Undefined function or method : nanewton" for input argumenTs of type» s=nane?rton< fun, dfun, 1» 1 亡亠5, io|o)(L 75040.73910. 73910.73910.7391if割线方式迭代法1. 程序：function x=ge_xian_fa(fname,dfname,x0,x1,e,N) %割线方式迭代法解

19、方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数， x0,x1 分别为迭代初值， e 为精度要求k=0;a=x1;b=x0;if nargin<5,N=500;endif nargin<4 e=1e-4;endwhile abs(a-b)>e&k<N,k=k+1;x=x1-(x1-x0)/(feval(fname,x1)-feval(fname,x0)*feval(fname,x1);if feval(fname,x)*feval(fname,x0)>0,x0=x;b=x0;elsex1=x;a=x1;endx=x1-(x1

20、-x0)/(feval(fname,x1)-feval(fname,x0)*feval(fname,x1);numb=int2str(k);disp('第 ',numb,' 次计算后 x=')fprintf('%fnn',x);endif k=N,warning（' 已达迭代次数上限 '）;end2. 运行结果：# Comma nd WindowFile Edit Debug Desktop Window Help_ X7 » fun=inline('x-cos(x)') ：dfun=inline(T+s

21、in(x)')：A5 » ge.xian.fa(fun, dfun, 1, 2, le-5, 100)=第1次计算后%=一e 0.742299第2次计算后:e 0.739491 1第3次计算后0. 739136第4次计算后乂=0. 739092第5次计算后-0.7390867X第6次计算后龙二 .0.739085L1第7次计算后沪 0.739085第8次计算后汶=0. 739085ii第9次计算后汶= 0.739085-第10次计算后瓦=二 0.739085V>4 Comma nd Window-XFile Edit Debug Desktop Window Helpu

22、. ZjyusoA第10次计算后:沪0. 7390854第11次计算后?c=0. 739085第12次计算后戸0. 739085第13次计算后沪0. 739085第14次计算后x-0. 739085第15次计算后:e0. 739085第16次计算后0. 739085第E次计算后:沪0. 739085ans =0. 7391»<>Newton 插值1. 程序：%保存文件名为 New_Int.m%Newton基本插值公式%x为向量，全部的插值节点%y为向量，差值节点处的函数值%xi 为标量，是自变量%yi 为 xi 出的函数估计值function yi=newton_chaz

23、hi(x,y,xi) n=length(x);m=length(y);if n=merror('The lengths of X ang Y must be equal!'); return;end%计算均差表 YY=zeros(n);Y(:,1)=y'for k=1:n-1for i=1:n-kif abs(x(i+k)-x(i)<eps error('the DATA is error!');return;endY(i,k+1)=(Y(i+1,k)-Y(i,k)/(x(i+k)-x(i);endend%计算牛顿插值公式yi=0;for i=1:

24、nz=1;for k=1:i-1z=z*(xi-x(k);endyi=yi+Y(1,i)*z;end2.运行结果:匚ommand WindowFile Edit Debug Desktop Window Help>> clear>> yi=newtoTi_chazhi ('pi/6 ; pi/4 t 0. 5 ; 0. 7071； f 2*pi/9)L? Undef mad function or method ' nevton_chashih for input arguments of ty -yi=newton_chazhi (:pi/S ; pi/41, 0.5:0. 70711,2*pi/9)0.6391yi=nMrttm_chazhi<pi/fi;pi/4:pi/3, 0 3：0 7071；0. S660：, 2*pi/9)OVRLagrange 插值1. 程序：function y0 = Language(x,y,x0) syms t l;if length(x)=length(y)n = length(x);elsedisp('x 和 y 的维数不相等！ ');return; % 检错end h=sym(0);for i=1:nl=sym(y(i);fo