数值分析编程及运行结果资料

1、高斯消元法1.程序: clear format ratA=input(' 输入增广矩阵 A=') m,n=size(A);for i=1:(m-1) numb=int2str(i);disp('第',numb,'次消元后的增广矩阵')for j=(i+1):m A(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:1 x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(

2、i,i); 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. 程序 :clearformat ratA=input(' 输入增广矩阵 A=') m,n=size(A);for i=1:(m-1)numb=int2str(i);disp('第'

3、;,n umb,'次选列主元后的增广矩阵') 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:1x(i)=(A

4、(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 DesktopWindow 旦e

5、lp0XXXVA第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='); f=inpu

6、t(' 输入增广矩阵最后一列元素 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 i=n-1:-

7、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(' 输入增广矩阵最后一列 b='

8、;); 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:n x(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 Mndow-XFile Edit Debug Desktop WindowHelp输入系数拒阵产4, 3,0;3,4,-1;-1, 0, 4A31 _43034-1-1

9、04输入增广矩阵最后一列2244山-24k =24ans =96/31120/31*152/31VOVR雅戈比迭代格式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:n y(i)=(b(i)-a(i,1:n)*x(1:n)+a(i,i)*x(i)/a(i,i); endfor i=1:nsum=sum+(y

10、(i)-x(i)八2;endif sqrt(sum)<ek break;elsefor i=1:n x(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（

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

12、;ww=input（' 输入松弛因子 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（' 迭代次

13、数为 '）;ndisp（' 方程组的解为 '）;xreturn;%上面：达到精度要求就结束程序，输出迭代次数和方程组 的解endX=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的 最终结果（并不是方程组的解）disp（'在最大迭代次数内不收敛！'）；disp（'最大迭代次数后的结果为'）; Command Window一XFil 亡 Edit D 亡 bug Desktop .Vindow Help»输入系数矩 HAt4. 3,0:3. 4.-1;-!. 0 4 输入方程组右端的列向Sb=-24

14、: 30 ;-24 输入迭代初值构戚的列向里U ：1：11： 输入最大迭代次数nni=10Ci输入误差精度尸山5e-输入松弛因子如戶1* 25迭代次敎为An =17方程组的解対X =3.0968|3. 8710-5*225BV<>OVR .：二分法求解方程的根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=x

15、elsea=xendk=k+1endf毅 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 ):

16、7; x=nabisect 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)%newton 迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数，x0为迭代初值，e为精度要求x=x0;x0=x+2*e;k=0;if nargin<

17、5,N=500;endif 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(' 已达迭代次数上限 ');end矶 Command WindowFile Edit Debug Desktop Jfindow HelpA, P"«S h>> fun-ml me ( xcos 'x.1' ) : dfunFinline ('

18、 1+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);end if k=N,warning（'已达迭代次数上限'）;end2.运行结果：匚ommand WindowFile Edit Debug Desktop Mndow Help» fun=inline x-cos i.x ) ： dfun=inlineC 1+s in (xj1 n): >&g

21、t; ge_xian_fa<fun, dfunT lj 2, le5t 100)第1次计算后0,742299第2次计算后艺二（L 739491第3次计算后XN01 739136第丄次计算后苴二0.739092第5次计算后0l 739096第占次计算后忙=0.739085第次计算后？0.739085第E次计算后忙=:0."39085IT第9次计算后艺二1 0.739035第10次计算后尸0i 7390854 Comma nd Window-XFile Edit Debug Desktop Window Helpu. ZjyusoA第10次计算后:沪0. 7390854第11次计算

22、后?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_chazhi(x,y,xi) n=length(x);m=

23、length(y);if n=merror('The lengths of X ang Y must be equal!'); return;end%计算均差表 Y Y=zeros(n);Y(:,1)=y'for k=1:n-1for i=1:n-k if 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:nz=1;for k=1:i-1z=z*(xi

24、-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

25、)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);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;e