王能超計(jì)算方法——算法設(shè)計(jì)及MATLAB實(shí)現(xiàn)課后代碼

1、第一章 插值方法1.1 Lagrange插值1.2 逐步插值1.3 分段三次Hermite插值1.4 分段三次樣條插值第二章 數(shù)值積分2.1 Simpson公式2.2 變步長(zhǎng)梯形法2.3 Romberg加速算法2.4 三點(diǎn)Gauss公式第三章 常微分方程德差分方法3.1 改進(jìn)的Euler方法3.2 四階Runge-Kutta方法3.3 二階Adams預(yù)報(bào)校正系統(tǒng)3.4 改進(jìn)的四階Adams預(yù)報(bào)校正系統(tǒng)第四章 方程求根4.1 二分法4.2 開(kāi)方法4.3 Newton下山法4.4 快速弦截法第五章 線性方程組的迭代法5.1 Jacobi迭代5.2 Gauss-Seidel迭代5.3 超松弛迭代5.

2、4 對(duì)稱超松弛迭代第六章 線性方程組的直接法6.1 追趕法6.2 Cholesky方法6.3 矩陣分解方法6.4 Gauss列主元消去法第一章 插值方法1.1 Lagrange插值計(jì)算Lagrange插值多項(xiàng)式在x=x0處的值.MATLAB文件：（文件名：Lagrange_eval.m） function y0,N= Lagrange_eval(X,Y,x0) %X,Y是已知插值點(diǎn)坐標(biāo) %x0是插值點(diǎn) %y0是Lagrange插值多項(xiàng)式在x0處的值 %N是Lagrange插值函數(shù)的權(quán)系數(shù) m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; fo

3、r j=1:m if j=i; N(i)=N(i)*(x0-X(j)/(X(i)-X(j); end end y0=y0+Y(i)*N(i);end用法X=;Y=;x0= ;y0,N= Lagrange_eval(X,Y,x0)1.2 逐步插值計(jì)算逐步插值多項(xiàng)式在x=x0處的值.MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0)%X,Y是已知插值點(diǎn)坐標(biāo)%x0是插值點(diǎn)%y0是Neville逐步插值多項(xiàng)式在x0處的值m=length(X);P=zeros(m,1);P1=zeros(m,1);P=Y;for i=1:m P1=

4、P; k=1; for j=i+1:m k=k+1; P(j)=P1(j-1)+(P1(j)-P1(j-1)*(x0-X(k-1)/(X(j)-X(k-1); endif abs(P(m)-P(m-1)<10-6; y0=P(m); return; endendy0=P(m);用法X=;Y=;x0= ;y0= Neville_eval(X,Y,x0)1.3 分段三次Hermite插值利用分段三次Hermite插值計(jì)算插值點(diǎn)處的函數(shù)近似值.MATLAB文件：（文件名：Hermite_interp.m）function y0=Hermite_interp(X,Y,DY,x0)%X,Y是已知插

5、值點(diǎn)向量序列%DY是插值點(diǎn)處的導(dǎo)數(shù)值%x0插值點(diǎn)橫坐標(biāo)%y0為待求的分段三次Hermite插值多項(xiàng)式在x0處的值%N向量長(zhǎng)度N=length(X);for i=1:N if x0>=X(i)&& x0<=X(i+1) k=i; break; endenda1=x0-X(k+1);a2=x0-X(k);a3= X(k)- X(k+1);y0=(a1/a3)2*(1-2*a2/a3)*Y(k)+(-a2/a3)2*(1+2*a1/a3)*Y(k+1)+. (a1/a3)2*a2*DY(k)+(-a2/a3)2*a1*DY(k+1);用法X=;Y=關(guān)于X的函數(shù);DY=Y;

6、x0=插值橫坐標(biāo);y0=Hermite_interp(X,Y,DY,x0)1.4 分段三次樣條插值計(jì)算在插值點(diǎn)處的函數(shù)值，并用來(lái)擬合曲線.MATLAB文件：（文件名：Spline_interp.m）function y0,C=Spline_interp(X,Y,s0,sN,x0)%X,Y是已知插值坐標(biāo)%s0,sN是兩端點(diǎn)的一次導(dǎo)數(shù)值%x0是插值點(diǎn)%y0三次樣條函數(shù)在x0處的值%C是分段三次樣條函數(shù)的系數(shù)N=length(X);C=zeros(4,N-1); h=zeros(1,N-1);mu=zeros(1,N-1); lmt=zeros(1,N-1);d=zeros(1,N); %d表示右端

7、函數(shù)值h=X(1,2:N)-X(1,1:N-1);mu(1,N-1)=1; lmt(1,1)=1;mu(1,1:N-2)=h(1,1:N-2)/(h(1,1:N-2)+h(1,2:N-1);lmt(1,2:N-1)=h(1,2:N-1)/(h(1,1:N-2)+h(1,2:N-1);d(1,1)=6*(Y(1,2)-Y(1,1)/h(1,1)-s0)/h(1,1);d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1)/h(1,N-1)/h(1,N-1);d(1,2:N-1)=6*(Y(1,3:N)-Y(1,2:N-1)/h(1,2:N-1)(Y(1,2:N-1)-Y(1,1:N-2)/h

8、(1,1:N-2)/(h(1,1:N-2)+h(1,2:N-1);%追趕法解三對(duì)角方程組bit=zeros(1,N-1);bit(1,1)=lmt(1,1)/2;for i=2:N-1 bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1);endy=zeros(1,N);y(1,1)=d(1,1)/2;for i=2:N y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1)/(2-mu(1,i-1)*bit(1,i-1);endx=zeros(1,N);x(1,N)=y(1,N);for i=N-1:-1:1 x(1,i)=y(1,i)-bit(1,

9、i)*x(1,i+1);endv=zeros(1,N-1);v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1)/h(1,1:N-1);C(4,:)=Y(1,1:N-1);C(3,:)=v-h*(2*x(1,1:N-1)+x(1,2:N)/6;C(2,:)=x(1,1:N-1)/2;C(1,:)=(x(1,2:N)-x(1,1:N-1)/(6*h);if nargin<5 y0=0;else for j=1:N-1 if x0>=X(1,j)& x0<X(1,j+1) omg=x0-X(1,j); y0=(C(4,j)*omg+C(3,j)omg+C(2,j

10、)*omg+C(1,j); end endend算例1：給定數(shù)據(jù)表：Xi0.250.300.390.450.53yi0.50000.54770.62450.67080.7280試求三次樣條插值函數(shù)S(x)，其中：S(0.25)=1.0000,S(0.53)=0.6868. 解：X=0.25,0.30,0.39,0.45,0.53;Y=0.5000,0.5477,0.6254,0.6708,0.7280;s0=1.0000;sN=0.6868;y0,C=Spline_interp(X,Y,s0,sN,x0);plot(0.25:0.01:0.30,polyval(C(:.1),0:0.01:0.

11、05),r-.);hold onplot(0.30:0.01:0.39,polyval(C(:.2),0:0.01:0.09),b);plot(0.39:0.01:0.45,polyval(C(:.3),0:0.01:0.06),k-*);plot(0.45:0.01:0.53,polyval(C(:.4),0:0.01:0.08);第二章 數(shù)值積分2.1 Simpson公式利用復(fù)化Simpson公式求被積函數(shù)f(x)在給定區(qū)間上的積分值.MATLAB文件：（文件名：FSimpson.m） function S=FSimpson(f,a,b,N) %f表示被積函數(shù)句柄 %a,b表示被積區(qū)間端點(diǎn)

12、 %N表示區(qū)間個(gè)數(shù) %S是用復(fù)化Simpson公式求得的積分值h=(b-a)/N;fa=feval(f,a);fb=feval(f,b);S=fa+fb;x=a;for i=1:N x=x+h/2; fx=feval(f,x); S=S+4*fx; x=x+h/2; fx=feval(f,x); S=S+2*fx;endS=h*S/6;算例2：利用復(fù)化Simpson公式計(jì)算積分.解：后面都要用到的f1:function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;運(yùn)行S=FSimpson(f,a,b,N)這里的N值需要自己輸入。2.2 變步長(zhǎng)梯形法利用變步長(zhǎng)梯形法求被積函數(shù)

13、f(x)在給定區(qū)間上的積分值.MATLAB文件：（文件名：bbct.m） function T,n=bbct(f,a,b,eps)%f表示被積函數(shù)句柄%a,b被積區(qū)間端點(diǎn)%eps精度%T是用變步長(zhǎng)梯形法求得的積分值%n表示二分區(qū)間的次數(shù)h=b-a;fa=feval(f,a);fb=feval(f,b);T1=h*(fa+fb)/2;T2=T1/2+h*feval(f,a+h/2)/2;n=1;%按變步長(zhǎng)梯形法求積分值while abs(T2-T1)>=eps h=h/2; T1=T2; S=0; x=a+h/2; while x<b fx=feval(f,x); S=S+fx; x

14、=x+h; end T2=T1/2+S*h/2; n=n+1;endT=T2;算例3：利用變步長(zhǎng)梯形法計(jì)算積分.解：function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;運(yùn)行 T,n=bbct(f,a,b,eps)這里的eps值需要自己輸入。2.3 Romberg加速算法利用Romberg加速算法計(jì)算被積函數(shù)f(x)在給定區(qū)間上的積分值.MATLAB文件：（文件名：Romberg.m） function quad,R=Romberg(f,a,b,eps)%f表示被積函數(shù)句柄%a,b被積區(qū)間端點(diǎn)%eps精度%quad是用Romberg加速算法求得的積分值%R為Romb

15、erg表%err表示誤差的估計(jì)h=b-a;R(1,1)=h*(feval(f,a)+feval(f,b)/2;M=1; J=0; err=1;while err>eps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1);endquad=R(J+1,J+1);算例4：利用Romber

16、g加速算法計(jì)算積分.解：function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;運(yùn)行 quad,R=Romberg(f,a,b,eps)這里的eps值需要自己輸入。2.4 三點(diǎn)Gauss公式利用三點(diǎn)Gauss公式計(jì)算被積函數(shù)f(x)在給定區(qū)間上的積分值.MATLAB文件：（文件名：TGauss.m） function G=TGauss(f,a,b)%f表示被積函數(shù)句柄%a,b被積區(qū)間端點(diǎn)%G是用三點(diǎn)Gauss公式求得的積分值x1=(a+b)/2-sqrt(3/5)*(b-a)/2;x2=(a+b)/2+sqrt(3/5)*(b-a)/2;G=(b-a)*(5*feva

17、l(f,x1)/9+8*feval(f,(a+b)/2)/9+5*feval(f,x2)/9)/2;算例5：利用三點(diǎn)Gauss公式計(jì)算積分.解：function f=f1(x)f=x/(4+x2);令f=f1;a=0;b=1;運(yùn)行 G=TGauss(f,a,b)第三章 常微分方程德差分方法3.1 改進(jìn)的Euler方法用改進(jìn)的Euler方法求解常微分方程.MATLAB文件：（文件名：MendEuler.m） function E=MendEuler(f,a,b,N,ya)%f是微分方程右端函數(shù)句柄%a,b是自變量的取值區(qū)間a,b的端點(diǎn)%N是區(qū)間等分的個(gè)數(shù)%ya表初值y(a)%E=x,y是自變量X

18、和解Y所組成的矩陣h=(b-a)/N;y=zeros(1,N+1);x=zeros(1,N+1);y(1)=ya;x=a:h:b;for i=1:N y1=y(i)+h*feval(f,x(i),y(i); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2;endE=x,y; 算例6：對(duì)于微分方程，用改進(jìn)的Euler方法求解.解：function z=f2(x,y)z=x2-y; 為衡量數(shù)值解的精度，我們求出該方程的解析解為.在此也以文件的形式表示如下：function y=solvef2(x)y=-exp(-x)+x2-2*x+2;令f=f2; a

19、=0; b=1; N=10; ya=1;運(yùn)行：E=MendEuler(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=E,y其中m內(nèi)的y表示準(zhǔn)確解.3.2 四階Runge-Kutta方法用四階Runge-Kutta方法求解常微分方程.MATLAB文件：（文件名：RungKutta4.m） function R=RungKutta4(f,a,b,N,ya)%f是微分方程右端函數(shù)句柄%a,b是自變量的取值區(qū)間a,b的端點(diǎn)%N是區(qū)間等分的個(gè)數(shù)%ya表初值y(a)%R=x,y是自變量X和解Y所組成的矩陣h=(b-a)/N;x=zeros(1,N+1);y=zeros(1,

20、N+1);x=a:h:b;y(1)=ya;for i=1:N k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4);endR=x,y; 算例7：對(duì)于微分方程，用四階Runge-Kutta方法求解.解：function z=f2(x,y)z=x2-y; 為衡量數(shù)值解的精度，我們求出該方程的解析解為.在此也以文件的形式表示如下：func

21、tion y=solvef2(x)y=-exp(-x)+x2-2*x+2;令f=f2; a=0; b=1; N=10; ya=1;運(yùn)行：R=RungKutta4(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=R,y其中m內(nèi)的y表示準(zhǔn)確解,精度比前者高.3.3 二階Adams預(yù)報(bào)校正系統(tǒng)用二階Adams預(yù)報(bào)校正系統(tǒng)求解常微分方程.MATLAB文件：（文件名：Adams2PC.m） function A=Adams2PC(f,a,b,N,ya)%f是微分方程右端函數(shù)句柄%a,b是自變量的取值區(qū)間a,b的端點(diǎn)%N是區(qū)間等分的個(gè)數(shù)%ya表初值y(a)%A=x,y是自變量

22、X和解Y所組成的矩陣h=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;for i=1:N if i=1 y1=y(i)+h*feval(f,x(i),y(i); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; dy1=feval(f,x(i),y(i); dy2=feval(f,x(i+1),y(i+1); else y(i+1)=y(i)+h*(3*dy2-dy1)/2; P=feval(f,x(i+1),y(i+1);l y(i+1)=y(i)+h*(P+dy2)/2; dy1=d

23、y2; dy2=feval(f,x(i+1),y(i+1); end endAx,y; 3.4 改進(jìn)的四階Adams預(yù)報(bào)校正系統(tǒng)用改進(jìn)的四階Adams預(yù)報(bào)校正系統(tǒng)求解常微分方程.MATLAB文件：（文件名：CAdams4PC.m） function A=CAdams4PC(f,a,b,N,ya)%f是微分方程右端函數(shù)句柄%a,b是自變量的取值區(qū)間a,b的端點(diǎn)%N是區(qū)間等分的個(gè)數(shù)%ya表初值y(a)%A=x,y是自變量X和解Y所組成的矩陣if N<4 break;endh=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;F=zer

24、os(1,4);for i=1:N if i<4 %用四階Runge-Kutta方法求初始解 k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); elseif i=4 F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); %預(yù)報(bào) p=feval(f,x(

25、i+1),py); F=F(2) F(3) F(4) p; y(i+1)=y(i)+(h/24)*(F*1,-5,19,9); %校正 p=py; c=y(i+1); else F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); %預(yù)報(bào) my=py-251*(p-c)/270; %改進(jìn) m=feval(f,x(i+1),my); F=F(2) F(3) F(4) m; cy=y(i)+(h/24)*(F*1,-5,19,9); %校正 y(i+1)=cy+19*(py-cy)/270; %改進(jìn) p=py; c=cy; en

26、dend A=x,y;附件(補(bǔ)充)：四階Adams預(yù)報(bào)校正系統(tǒng)的程序：MATLAB文件：（文件名：Adams4PC.m） function A=Adams4PC(f,a,b,N,ya)%f是微分方程右端函數(shù)句柄%a,b是自變量的取值區(qū)間a,b的端點(diǎn)%N是區(qū)間等分的個(gè)數(shù)%ya表初值y(a)%A=x,y是自變量X和解Y所組成的矩陣if N<4 break;endh=(b-a)/N;x=zeros(1,N+1);y=zeros(1,N+1);x=a:h:b;y(1)=ya;F=zeros(1,4);for i=1:N if i<4 %用四階Runge-Kutta方法求初始解 k1=fev

27、al(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); else F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); %預(yù)報(bào) p=feval(f,x(i+1),py); F=F(2) F(3) F(4) p; y(i+1)=y(i)+(h/24)*(F*1,-5,19,9)

28、; %校正 endend A=x,y;算例8：分別用二階Adams預(yù)報(bào)校正系統(tǒng)，四階Adams預(yù)報(bào)校正系統(tǒng)和改進(jìn)四階Adams預(yù)報(bào)校正系統(tǒng)求解如下微分方程初值問(wèn)題：.解：function z=f3(x,y)z=-y+x+1; 為衡量數(shù)值解的精度，我們求出該方程的解析解為.在此也以文件的形式表示如下：function y=solvef3(x)y=exp(-x)+x;令f=f3; a=0; b=1; N=10; ya=1;運(yùn)行：A2=Adams2PC(f,a,b,N,ya)； A4=Adams4PC(f,a,b,N,ya); CA4=CAdams4PC(f,a,b,N,ya); y=solvef3

29、(a:(b-a)/N:b); m=A2,A4(:,2),CA4(:.2),y其中m共有5列數(shù)據(jù)：左右依次為離散節(jié)點(diǎn)值、二階Adams預(yù)報(bào)校正系統(tǒng)所求解、四階Adams預(yù)報(bào)校正系統(tǒng)所求解、改進(jìn)四階Adams預(yù)報(bào)校正系統(tǒng)所求解和準(zhǔn)確解,精度依次遞增.第四章 方程求根4.1 二分法用二分法求解非線性方程在區(qū)間a,b內(nèi)的根.MATLAB文件：（文件名：demimethod.m） function x,k=demimethod(a,b,f,emg)%a,b表示求解區(qū)間a,b的端點(diǎn)%f表示所求解方程德函數(shù)名%emg是精度指標(biāo)%x表示所求近似解%k表示循環(huán)次數(shù)fa=feval(f,a);fab=feval(

30、f,(a+b)/2);k=0;while abs(b-a)>emg if fab=0 x=(a+b)/2； return; elseif fa*fab<0 b=(a+b)/2; else a=(a+b)/2; end fa=feval(f,a); fab=feval(f,(a+b)/2); k=k+1;endx=(a+b)/2;算例9：求解方程在區(qū)間0,內(nèi)的實(shí)根，使精度達(dá)到.解：function f=func2(x) f=sqrt(x2+1)-tan(x); 輸入：f=func2; a=0; b=pi/2; emg=10-5; x0,k=demimethod(a,b,f,emg)4

31、.2 開(kāi)方法求實(shí)數(shù)的開(kāi)方運(yùn)算.MATLAB文件：（文件名：Kaifang.m） function y=Kaifang(a,eps,x0)%a是被開(kāi)方數(shù)%eps是精度指標(biāo)%x0表示初值%y是a的開(kāi)方x(1)=x0;x(2)=(x(1)+a/x(1)/2;k=2;while abs(x(k)-x(k-1)>eps x(k+1)=(x(k)+a/x(k)/2; k=k+1;endy=x;算例10：用開(kāi)方法求，設(shè)取x0=1.解：令a=2; eps=10-6; x0=1; 運(yùn)行 y=Kaifang(a,eps,x0)4.3 Newton下山法用牛頓下山法求解非線性方程的根.MATLAB文件：（文件

32、名：Mendnewton.m） function x,k=Mendnewton(f,x0,emg)%用牛頓下山法求解非線性方程%f表示非線性方程%x0是迭代初值，此方法是局部收斂，初值選擇要恰當(dāng)%emg是精度指標(biāo)%k,u分別表示迭代次數(shù)和下山因子f1,d1=feval(f,x0); %d1表示非線性方程f在x0處的導(dǎo)數(shù)值，以下類同k=1;x(1)=x0;x(2)=x(1)-f1/d1;while abs(f1)>emg %控制精度 u=1; k=k+1; f1,d1=feval(f,x(k); x(k+1)=x(k)-u*f1/d1; %牛頓下山迭代 f2,d2=feval(f,x(k+

33、1); while abs(f2)>abs(f1) %保證迭代后的函數(shù)值比迭代前的函數(shù)值小 u=u/2; x(k+1)=x(k)-u*f1/d1; %牛頓下山迭代 f2,d2=feval(f,x(k+1); endend算例11:用牛頓下山法求方程的根，使精度達(dá)到.初值分別為：(1) x0=-1.2; (2)x0=2.0;解：編寫func3.m.function f,d=func3(x)f=sqrt(x2+1)-tan(x);d1=sqrt(x2+1)-tan(x);d=subs(diff(d1); %對(duì)函數(shù)f求一次導(dǎo)數(shù)在命令窗口中輸入： f=func3; x,k=Mendnewton(

34、f,x0,10-6) 其中，x0需要自己輸入.4.4 快速弦截法用快速弦截法求解非線性方程的根.MATLAB文件：（文件名：Fast_chord.m） function x,k=Fast_chord(f,x1,x2,emg)%用快速弦截法求解非線性方程%f表示非線性方程函數(shù)%x1,x2是迭代初值%emg是精度指標(biāo)%k表示循環(huán)次數(shù)k=1;y1=feval(f,x1);y2=feval(f,x2);x(k)=x2-(x2-x1)*y2/(y2-y1); %用快速弦截法進(jìn)行迭代求解y(k)=feval(f,x(k);k=k+1;x(k)=x(k-1)-(x(k-1)-x2)*y(k-1)/(y(k-

35、1)-y2);while abs(x(k)-x(k-1)>emg %控制精度 y(k)=feval(f,x(k); x(k+1)=x(k)-(x(k)-x(k-1)*y(k)/(y(k)-y(k-1); k=k+1;end算例12:用快速弦截法求解非線性方程的根，要求精度為，初值為：解：編寫函數(shù)文件func4.m function f=func4(x) f=exp(x)-4*cos(x)在命令窗口輸入： f=func4; x,k=Fast_chord(f,pi/4,pi/2,10-6)第五章 線性方程組的迭代法5.1 Jacobi迭代用Jacobi迭代法求解線性方程組.MATLAB文件：

36、（文件名：Jacobimethod.m） function x,k=Jacobimethod(A,b,x0,N,emg)%A是線性方程組的左端矩陣%b是右端向量%x0是迭代初始值%N表示迭代次數(shù)上限，若迭代次數(shù)大于N，則迭代失敗%emg表示控制精度%用Jacobi迭代法求線性方程組A*x=b的解%k表示迭代次數(shù)%x表示用迭代法求得的線性方程組的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while r>emg for i=1:n sum=0; for j=1:n if i=j sum=s

37、um+A(i,j)*x1(j); end end x2(i)=(b(i)-sum)/A(i,i); end r=max(abs(x2-x1); x1=x2; k=K+1; if k>N disp(迭代失敗，返回); return; endendx=x1;算例13：用Jacobi迭代法求解方程組： 其精確解為x=-1,-1,-1,-1.解：令A(yù)=-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4; b=1,1,1,1; x0=0,0,0,0;在命令窗口中輸入： x,k=Jacobimethod(A,b,x0,100,10-5)5.2 Gauss-Seidel迭代用Gau

38、ss-Seidel迭代法求解線性方程組.MATLAB文件：（文件名：Gaussmethod.m） function x,k=Gaussmethod(A,b,x0,N,emg)%A是線性方程組的左端矩陣%b是右端向量%x0是迭代初始值%N表示迭代次數(shù)上限，若迭代次數(shù)大于N，則迭代失敗%emg表示控制精度%用Gauss-Seidel迭代法求線性方程組A*x=b的解%k表示迭代次數(shù)%x表示用迭代法求得的線性方程組的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while r>emg for i=

39、1:n sum=0; for j=1:n if j>i sum=sum+A(i,j)*x1(j); elseif j<i sum=sum+A(i,j)*x2(j); end end x2(i)=(b(i)-sum)/A(i,i); end r=max(abs(x2-x1); x1=x2; k=K+1; if k>N disp(迭代失敗，返回); return; endendx=x1;算例14：用Gauss-Seidel迭代法求解方程組： 其精確解為x=-1,-1,-1,-1.解：令A(yù)=-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4; b=1,1,1,1

40、; x0=0,0,0,0;在命令窗口中輸入： x,k=Gaussmethod(A,b,x0,100,10-5)5.3 超松弛迭代用超松弛(SOR)迭代法求解線性方程組.MATLAB文件：（文件名：SORmethod.m） function x,k=SORmethod(A,b,x0,N,emg,w)%A是線性方程組的左端矩陣%b是右端向量%x0是迭代初始值%N表示迭代次數(shù)上限，若迭代次數(shù)大于N，則迭代失敗%emg表示控制精度%w表示松弛因子%用SOR迭代法求線性方程組A*x=b的解%k表示迭代次數(shù)%x表示用迭代法求得的線性方程組的近似解n=length(A);x1=zeros(n,1); x2=

41、zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while r>emg for i=1:n sum=0; for j=1:n if j>=i sum=sum+A(i,j)*x1(j); elseif j<i sum=sum+A(i,j)*x2(j); end end x2(i)=x1(i)+w*(b(i)-sum)/A(i,i); end r=max(abs(x2-x1); x1=x2; k=K+1; if k>N disp(迭代失敗，返回); return; endendx=x1;算例15：用超松弛(SOR)迭代法求解方程組： 其精確解

42、為x=-1,-1,-1,-1.解：令A(yù)=-4,1,1,1;1,-4,1,1;1,1,-4,1;1,1,1,-4; b=1,1,1,1; x0=0,0,0,0;在命令窗口中輸入： x,k=SORmethod(A,b,x0,100,10-5,1)當(dāng)松弛因子為1時(shí)，超松弛迭代法等同于Gauss-Seidel迭代法.5.4 對(duì)稱超松弛迭代用對(duì)稱超松弛(SSOR)迭代法求解線性方程組.MATLAB文件：（文件名：SSORmethod.m） function x,k=SSORmethod(A,b,x0,N,emg,w)%A是線性方程組的左端矩陣%b是右端向量%x0是迭代初始值%N表示迭代次數(shù)上限，若迭代次

43、數(shù)大于N，則迭代失敗%emg表示控制精度%w表示松弛因子%用SSOR迭代法求線性方程組A*x=b的解%k表示迭代次數(shù)%x表示用迭代法求得的線性方程組的近似解n=length(A);x1=zeros(n,1); x2=zeros(n,1); x3=zeros(n,1);x1=x0; r=max(abs(b-A*x1);k=0;while r>emg for i=1:n sum=0; for j=1:n if j>i sum=sum+A(i,j)*x1(j); elseif j<i sum=sum+A(i,j)*x2(j); end end x2(i)=(1-w)*x1(i)+w

44、*(b(i)-sum)/A(i,i); end for i=n:-1:1 sum=0; for j=1:n if j>i sum=sum+A(i,j)*x3(j); elseif j<i sum=sum+A(i,j)*x2(j); end end x3(i)=(1-w)*x2(i)+w*(b(i)-sum)/A(i,i); end r=max(abs(x3-x1); x1=x3; k=K+1; if k>N disp(迭代失敗，返回); return; endendx=x1;算例16：用對(duì)稱超松弛(SSOR)迭代法和超松弛(SSOR)迭代法求解線性方程組： 精度控制在.解：令

45、A=4,-1,0;-1,4,-1;0,-1,4; b=1,4,-3; x0=0,0,0;在命令窗口中輸入： x1,k1=SORmethod(A,b,x0,100,10-5,1.2) x2,k2=SSORmethod(A,b,x0,100,10-5,1.2)第六章 線性方程組的直接法6.1 追趕法用追趕法解三對(duì)角線性方程組.MATLAB文件：（threedia.m） function x=threedia(a,b,c,f)%求解線性方程組Ax=f，其中A是三對(duì)角陣%a是矩陣A的下對(duì)角線元素a(1)=0%b是矩陣A的對(duì)角線元素%c是矩陣A的上對(duì)角線元素c(N)=0%f是方程組的右端向量N=leng

46、th(f);x=zeros(1,N); y=zeros(1,N); d=zeros(1,N); u=zeros(1,N);%預(yù)處理d(1)=b(1);for i=1:N-1 u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i+1)*u(i);end %追的過(guò)程y(1)=f(1)/d(1);for i=2:N y(i)=(f(i)-a(i)*y(i-1)/d(i);end%趕的過(guò)程x(N)=y(N);for i=N-1:-1:1 x(i)=y(i)-u(i)*x(i+1);end算例17：用追趕法求解方程組： 解：令a=0,-1,-1,-3; b=2,3,2,5; c=-1,-2

47、,-1,0; f=6,1,0,1;在命令窗口輸入： x=threedia(a,b,c,f)6.2 Cholesky方法用Cholesky分解法解對(duì)稱方程組.MATLAB文件：（文件名：Chol_decompose.m） function x=Chol_decompose(A,b)%用Cholesky分解求解線性方程組Ax=b%A是對(duì)稱矩陣%L是單位下三角矩陣%D是對(duì)角陣%對(duì)矩陣A進(jìn)行三角分解:A=LDLN=length(A);L=zeros(N,N); D=zeros(1,N);for i=1:N L(i,i)=1;endD(1)=A(1,1);for i=2:N for j=1:i-1 if j=1 L(i,j)=A(i,j)/D(j); else sum1=0; for k=1