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

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.運行結(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. 程序 :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.運行結(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 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.運行結(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(' 輸入增廣矩陣最后一

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.運行結(jié)果: Command Window 口File Edit Debug Desktop Window Help輸入系數(shù)拒陣尸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(' 輸入系數(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: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.運行結(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(A,b,X,nm,

11、w,ww)%用超松弛迭代法求解方程組 Ax=b %俞入：A為方程組的系數(shù)矩陣，b為方程組右端的列向量，X為迭 代初值構(gòu)成的列向量，nm為最大迭代次數(shù)，w為誤差精度，ww為松 弛因子%俞出： x 為求得的方程組的解構(gòu)成的列向量， n 為迭代次數(shù)A=input（' 俞入系數(shù)矩陣 A='）;b=input（' 俞入方程組右端的列向量 b='）;X=input（' 俞入迭代初值構(gòu)成的列向量 X='）;nm=input（' 俞入最大迭代次數(shù) 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;%計算迭代格式中的常數(shù)項%下面是迭代過程while n<=nmx=M*X+g; % 用迭代格式進行迭代if norm（x-X,'inf'）<wdisp（' 迭代次數(shù)為 '）;nd

13、isp（' 方程組的解為 '）;xreturn;% 上面：達到精度要求就結(jié)束程序， 輸出迭代次數(shù)和方程組的 解endX=x;n=n+1;end%下面：如果達到最大迭代次數(shù)仍不收斂， 輸出警告語句及迭代的最 終結(jié)果（并不是方程組的解）disp（' 在最大迭代次數(shù)內(nèi)不收斂 !'）;disp（' 最大迭代次數(shù)后的結(jié)果為 '）;2. 運行結(jié)果：d' Command WindowFile Edit Debug Desktop Window 旦elp »輸入系數(shù)矩P$A=：4, 3, 0:3, 4,-1 ;-l, 0, 4 輸入方程組右端的

14、列向堅=24： 30:-24 輸入迭代初值構(gòu)成的列向1X=1 ;1：1 輸入最大迭代次數(shù)nm=100輸入誤差精度w=0. 5e- 輸入松弛因子ww=l25 迭代次數(shù)為方程組的解為3. 09683.8710-5.2258v二分法求解方程的根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=xelsea=xendk=k+1 end2

15、.運行結(jié)果: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)及其導函數(shù)的M函數(shù)句柄或內(nèi)嵌函數(shù)，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(' 已達迭代次數(shù)上限 ');end2.運行結(jié)果:磯 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)及其導函數(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);endif k=N,warning（' 已達迭代次數(shù)上限 '）;end2. 運行結(jié)果：# 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為向量，全部的插值節(jié)點%y為向量，差值節(jié)點處的函數(shù)值%xi 為標量，是自變量%yi 為 xi 出的函數(shù)估計值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%計算牛頓插值公式y(tǒng)i=0;for i=1:

24、nz=1;for k=1:i-1z=z*(xi-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)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);fo