數(shù)值分析編程及運行結(jié)果資料

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運行結(jié)果:協(xié)匚ommand WindowWindowHelp-口XFile Edit DebugDesktop輸入増廣矩P5A=： 1116-13142-61 -5A -1116-13142-61-5第1次消元后的増廣矩陣A -1116042100-S-I'17第2次消元后的增廣矩時A =1116042100033回代求網(wǎng)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.運行結(jié)果：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 進(jìn)行分解 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.運行結(jié)果：和匚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(' 輸入系數(shù)矩陣 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.運行結(jié)果:機 Command Mndow-XFile Edit Debug Desktop WindowHelp輸入系數(shù)拒陣產(chǎn)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(' 輸入系數(shù)矩陣 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.運行結(jié)果:林匚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為方程組的系數(shù)矩陣，b為方程組右端的列向量， 代初值構(gòu)成的列向量， nm 為最大迭代次數(shù)， w 為誤差精度， 松弛因子%輸出： x 為求得的方程組的解構(gòu)成的列向量， n 為迭代次數(shù)A=input（' 輸入系數(shù)矩陣 A='）;b=i nput（'輸入方程組右端的列向量b='）;X=input（' 輸入迭代初值構(gòu)成的列向量 X='）;nm=input（' 輸入最大迭代次數(shù) 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;%計算迭代矩陣%計算迭代格式中的常數(shù)%令 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;% 用迭代格式進(jìn)行迭代if norm（x-X,'inf'）<wdisp（' 迭代次

13、數(shù)為 '）;ndisp（' 方程組的解為 '）;xreturn;%上面：達(dá)到精度要求就結(jié)束程序，輸出迭代次數(shù)和方程組 的解endX=x;n=n+1;end%下面：如果達(dá)到最大迭代次數(shù)仍不收斂，輸出警告語句及迭代的 最終結(jié)果（并不是方程組的解）disp（'在最大迭代次數(shù)內(nèi)不收斂！'）；disp（'最大迭代次數(shù)后的結(jié)果為'）; Command Window一XFil 亡 Edit D 亡 bug Desktop .Vindow Help»輸入系數(shù)矩 HAt4. 3,0:3. 4.-1;-!. 0 4 輸入方程組右端的列向Sb=-24

14、: 30 ;-24 輸入迭代初值構(gòu)戚的列向里U ：1：11： 輸入最大迭代次數(shù)nni=10Ci輸入誤差精度尸山5e-輸入松弛因子如戶1* 25迭代次敎為An =17方程組的解対X =3.0968|3. 8710-5*225BV<>OVR .：二分法求解方程的根1.程序：%其中 a， b 表示查找根存在的范圍， M 表示要求解函數(shù)的值 %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)及其導(dǎo)函數(shù)的M函數(shù)句柄或內(nèi)嵌函數(shù)，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(' 已達(dá)迭代次數(shù)上限 ');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)及其導(dǎo)函數(shù)的M函數(shù)句柄或內(nèi)嵌函 數(shù)，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（'已達(dá)迭代次數(shù)上限'）;end2.運行結(jié)果：匚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 為向量，全部的插值節(jié)點 %y 為向量，差值節(jié)點處的函數(shù)值 %xi 為標(biāo)量，是自變量 %yi 為 xi 出的函數(shù)估計值 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%計算牛頓插值公式y(tǒng)i=0;for i=1:nz=1;for k=1:i-1z=z*(xi

24、-x(k);endyi=yi+Y(1,i)*z;end2運行結(jié)果:匚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 的維數(shù)不相等！ ');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