王能超 计算方法 - 算法设计及MATLAB实现课后代码

第一章 插值方法 1.1 Lagrange插值 1.2 逐步插值

1.3 分段三次Hermite插值 1.4 分段三次样条插值

第二章 数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式

第三章 常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统

3.4 改进的四阶Adams预报校正系统

第四章 方程求根 4.1 二分法 4.2 开方法

4.3 Newton下山法 4.4 快速弦截法

第五章 线性方程组的迭代法 5.1 Jacobi迭代

5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代

第六章 线性方程组的直接法 6.1 追赶法

6.2 Cholesky方法 6.3 矩阵分解方法

6.4 Gauss列主元消去法

1

第一章 插值方法 1.1 Lagrange插值

计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m） function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点

%y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for 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 逐步插值

计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）

function y0=Neville_eval(X,Y,x0)

%X,Y是已知插值点坐标 %x0是插值点

%y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1;

for j=i+1:m k=k+1;

2

P(j)=P1(j-1)+(P1(j)-P1(j-1))*(x0-X(k-1))/(X(j)-X(k-1)); end

if abs(P(m)-P(m-1))<10^-6; y0=P(m); return; end end

y0=P(m);

用法》X=[…];Y=[…];

》x0= ;

》y0= Neville_eval(X,Y,x0)

1.3 分段三次Hermite插值

利用分段三次Hermite插值计算插值点处的函数近似值. MATLAB文件：（文件名：Hermite_interp.m）

function y0=Hermite_interp(X,Y,DY,x0) %X,Y是已知插值点向量序列 %DY是插值点处的导数值 %x0插值点横坐标

%y0为待求的分段三次Hermite插值多项式在x0处的值 %N向量长度 N=length(X); for i=1:N

if x0>=X(i)&& x0<=X(i+1) k=i; break; end end

a1=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=关于X的函数;DY=Y’;

》x0=插值横坐标;

》y0==Hermite_interp(X,Y,DY,x0)

1.4 分段三次样条插值

计算在插值点处的函数值，并用来拟合曲线. MATLAB文件：（文件名：Spline_interp.m）

function [y0,C]=Spline_interp(X,Y,s0,sN,x0)

3

%X,Y是已知插值坐标

%s0,sN是两端点的一次导数值 %x0是插值点

%y0三次样条函数在x0处的值 %C是分段三次样条函数的系数 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表示右端函数值 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(1,1:N-2))/(h(1,1:N-2)+h(1,2:N-1)); %追赶法解三对角方程组 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)); end

y=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)); end

x=zeros(1,N); x(1,N)=y(1,N); for i=N-1:-1:1

x(1,i)=y(1,i)-bit(1,i)*x(1,i+1); end

v=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

4

if x0>=X(1,j)& x0

y0=((C(4,j)*omg+C(3,j))omg+C(2,j))*omg+C(1,j); end end end

算例1：给定数据表： Xi 0.25 yi 0.5000 0.30 0.5477 0.39 0.6245 0.45 0.6708 0.53 0.7280 试求三次样条插值函数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.05),’r-.’); hold on

plot(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));

第二章 数值积分 2.1 Simpson公式

利用复化Simpson公式求被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：FSimpson.m） function S=FSimpson(f,a,b,N) %f表示被积函数句柄 %a,b表示被积区间端点 %N表示区间个数

%S是用复化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; end S=h*S/6;

5

算例2：利用复化Simpson公式计算积分S?解：后面都要用到的f1: function f=f1(x) f=x/(4+x^2);

令f=@f1;a=0;b=1;

运行S=FSimpson(f,a,b,N) 这里的N值需要自己输入。

xdx?04?x2.

12.2 变步长梯形法

利用变步长梯形法求被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：bbct.m） function [T,n]=bbct(f,a,b,eps) %f表示被积函数句柄 %a,b被积区间端点 %eps精度

%T是用变步长梯形法求得的积分值 %n表示二分区间的次数 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;

%按变步长梯形法求积分值 while abs(T2-T1)>=eps h=h/2; T1=T2; S=0; x=a+h/2; while x

fx=feval(f,x); S=S+fx; x=x+h; end

T2=T1/2+S*h/2; n=n+1; end T=T2;

算例3：利用变步长梯形法计算积分T?解：

6

xdx?04?x2.

1function f=f1(x) f=x/(4+x^2);

令f=@f1;a=0;b=1;

运行 [T,n]=bbct(f,a,b,eps)

这里的eps值需要自己输入。

2.3 Romberg加速算法

利用Romberg加速算法计算被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：Romberg.m） function [quad,R]=Romberg(f,a,b,eps) %f表示被积函数句柄 %a,b被积区间端点 %eps精度

%quad是用Romberg加速算法求得的积分值 %R为Romberg表 %err表示误差的估计 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))/(4^k-1); end

err=abs(R(J+1,J)-R(J+1,J+1)); end

quad=R(J+1,J+1);

算例4：利用Romberg加速算法计算积分R?解：

function f=f1(x) f=x/(4+x^2);

令f=@f1;a=0;b=1;

运行 [quad,R]=Romberg(f,a,b,eps) 这里的eps值需要自己输入。

7

xdx?04?x2.

1

2.4 三点Gauss公式

利用三点Gauss公式计算被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：TGauss.m） function G=TGauss(f,a,b) %f表示被积函数句柄 %a,b被积区间端点

%G是用三点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*feval(f,x1)/9+8*feval(f,(a+b)/2)/9+5*feval(f,x2)/9)/2;

算例5：利用三点Gauss公式计算积分R?解：

function f=f1(x) f=x/(4+x^2);

令f=@f1;a=0;b=1; 运行 G=TGauss(f,a,b)

第三章 常微分方程德差分方法 3.1 改进的Euler方法

用改进的Euler方法求解常微分方程. MATLAB文件：（文件名：MendEuler.m） function E=MendEuler(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%E=[x’,y’]是自变量X和解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; end

E=[x’,y’];

8

xdx?04?x2.

1

算例6：对于微分方程解：

function z=f2(x,y) z=x^2-y;

为衡量数值解的精度，我们求出该方程的解析解为y??e?x?x2?2x?2.在此也以文件的形式表示如下： function y=solvef2(x) y=-exp(-x)+x^2-2*x+2;

令f=@f2; a=0; b=1; N=10; ya=1; 运行：E=MendEuler(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=[E,y’]

其中m内的y’表示准确解.

3.2 四阶Runge-Kutta方法

用四阶Runge-Kutta方法求解常微分方程. MATLAB文件：（文件名：RungKutta4.m） function R=RungKutta4(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%R=[x’,y’]是自变量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

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); end

R=[x’,y’];

算例7：对于微分方程解：

9

dy?x2?y,y(0)?1,0?x?1，用改进的Euler方法求解. dxdy?x2?y,y(0)?1,0?x?1，用四阶Runge-Kutta方法求解. dxfunction z=f2(x,y) z=x^2-y;

为衡量数值解的精度，我们求出该方程的解析解为y??e?x?x2?2x?2.在此也以文件的形式表示如下： function y=solvef2(x) y=-exp(-x)+x^2-2*x+2;

令f=@f2; a=0; b=1; N=10; ya=1; 运行：R=RungKutta4(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=[R,y’]

其中m内的y’表示准确解,精度比前者高.

3.3 二阶Adams预报校正系统

用二阶Adams预报校正系统求解常微分方程. MATLAB文件：（文件名：Adams2PC.m） function A=Adams2PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量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=dy2;

dy2=feval(f,x(i+1),y(i+1)); end end A[x’,y’];

10

3.4 改进的四阶Adams预报校正系统

用改进的四阶Adams预报校正系统求解常微分方程. MATLAB文件：（文件名：CAdams4PC.m） function A=CAdams4PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量X和解Y所组成的矩阵 if N<4 break; end

h=(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=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]’); %预报 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]’); %校正 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]’); %预报 my=py-251*(p-c)/270; %改进 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; %改进 p=py; c=cy; end

11

end A=[x’,y’];

附件(补充)：四阶Adams预报校正系统的程序： MATLAB文件：（文件名：Adams4PC.m） function A=Adams4PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量X和解Y所组成的矩阵 if N<4 break; end

h=(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=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); else

F=feval(f,x(i-3:i),y(i-3:i));

py=y(i)+(h/24)*(F*[-9,37,-59,55]’); %预报 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]’); %校正 end end A=[x’,y’];

算例8：分别用二阶Adams预报校正系统，四阶Adams预报校正系统和改进四阶Adams预报校正系统求解如下微分方程初值问题：解：

function z=f3(x,y) z=-y+x+1;

12

dy??y?x?1,y(0)?1,0?x?1. dx 为衡量数值解的精度，我们求出该方程的解析解为y?e?x?x.在此也以文件的形式表示如下：

function y=solvef3(x) y=exp(-x)+x;

令f=@f3; a=0; b=1; N=10; ya=1; 运行：A2=Adams2PC(f,a,b,N,ya)； A4=Adams4PC(f,a,b,N,ya); CA4=CAdams4PC(f,a,b,N,ya); y=solvef3(a:(b-a)/N:b); m=[A2,A4(:,2),CA4(:.2),y’]

其中m共有5列数据：左→右依次为离散节点值、二阶Adams预报校正系统所求解、四阶Adams预报校正系统所求解、改进四阶Adams预报校正系统所求解和准确解,精度依次递增.

第四章 方程求根 4.1 二分法

用二分法求解非线性方程f(x)?0在区间[a,b]内的根. MATLAB文件：（文件名：demimethod.m） function [x,k]=demimethod(a,b,f,emg) %a,b表示求解区间[a,b]的端点 %f表示所求解方程德函数名 %emg是精度指标 %x表示所求近似解 %k表示循环次数 fa=feval(f,a);

fab=feval(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; end

x=(a+b)/2;

13

算例9：求解方程f(x)?解：

function f=func2(x)

f=sqrt(x^2+1)-tan(x);

x2?1?tanx在区间[0,?/2]内的实根，使精度达到10?5.

输入：f=@@func2; a=0; b=pi/2; emg=10^-5; [x0,k]=demimethod(a,b,f,emg)

4.2 开方法

求实数的开方运算. MATLAB文件：（文件名：Kaifang.m） function y=Kaifang(a,eps,x0) %a是被开方数 %eps是精度指标 %x0表示初值 %y是a的开方 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; end y=x’;

算例10：用开方法求2，设取x0=1. 解：令a=2; eps=10^-6; x0=1; 运行 y=Kaifang(a,eps,x0)

4.3 Newton下山法

用牛顿下山法求解非线性方程f(x)?0的根.

MATLAB文件：（文件名：Mendnewton.m） function [x,k]=Mendnewton(f,x0,emg) %用牛顿下山法求解非线性方程 %f表示非线性方程

%x0是迭代初值，此方法是局部收敛，初值选择要恰当 %emg是精度指标

%k,u分别表示迭代次数和下山因子

[f1,d1]=feval(f,x0); ?表示非线性方程f在x0处的导数值，以下类同 k=1; x(1)=x0;

14

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+1));

while abs(f2)>abs(f1) %保证迭代后的函数值比迭代前的函数值小 u=u/2;

x(k+1)=x(k)-u*f1/d1; %牛顿下山迭代 [f2,d2]=feval(f,x(k+1)); end end

算例11:用牛顿下山法求方程f(x)?(1)x0=-1.2; (2)x0=2.0;

解：编写func3.m. function [f,d]=func3(x) f=sqrt(x^2+1)-tan(x); d1=’sqrt(x^2+1)-tan(x)’;

d=subs(diff(d1)); %对函数f求一次导数 在命令窗口中输入：

f=@func3; [x,k]=Mendnewton(f,x0,10^-6) 其中，x0需要自己输入.

4.4 快速弦截法

用快速弦截法求解非线性方程f(x)?0的根. MATLAB文件：（文件名：Fast_chord.m） function [x,k]=Fast_chord(f,x1,x2,emg)

%用快速弦截法求解非线性方程 %f表示非线性方程函数 %x1,x2是迭代初值 %emg是精度指标 %k表示循环次数 k=1;

y1=feval(f,x1); y2=feval(f,x2);

x(k)=x2-(x2-x1)*y2/(y2-y1); %用快速弦截法进行迭代求解 y(k)=feval(f,x(k)); k=k+1;

x(k)=x(k-1)-(x(k-1)-x2)*y(k-1)/(y(k-1)-y2);

15

x2?1?tanx的根，使精度达到10?6.初值分别为：

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:用快速弦截法求解非线性方程4cosx?e的根，要求精度为??10，初值为：

x?6x1??/4,x2??/2.

解：编写函数文件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文件：（文件名：Jacobimethod.m） function [x,k]=Jacobimethod(A,b,x0,N,emg)

%A是线性方程组的左端矩阵 %b是右端向量 %x0是迭代初始值

%N表示迭代次数上限，若迭代次数大于N，则迭代失败 %emg表示控制精度

%用Jacobi迭代法求线性方程组A*x=b的解 %k表示迭代次数

%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=sum+A(i,j)*x1(j); end end

x2(i)=(b(i)-sum)/A(i,i);

16

end

r=max(abs(x2-x1)); x1=x2; k=K+1; if k>N

disp(‘迭代失败，返回’); return; end end x=x1;

算例13：用Jacobi迭代法求解方程组：

??4x1?x2?x3?x4?1?x?4x?x?x?1?1234 ?

?x1?x2?4x3?x4?1??x1?x2?x3?4x4?1其精确解为x=[-1,-1,-1,-1]’.

解：令A=[-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迭代

用Gauss-Seidel迭代法求解线性方程组. MATLAB文件：（文件名：Gaussmethod.m） function [x,k]=Gaussmethod(A,b,x0,N,emg)

%A是线性方程组的左端矩阵 %b是右端向量 %x0是迭代初始值

%N表示迭代次数上限，若迭代次数大于N，则迭代失败 %emg表示控制精度

%用Gauss-Seidel迭代法求线性方程组A*x=b的解 %k表示迭代次数

%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

17

if j>i

sum=sum+A(i,j)*x1(j); elseif j

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; end end x=x1;

算例14：用Gauss-Seidel迭代法求解方程组：

??4x1?x2?x3?x4?1?x?4x?x?x?1?1234 ?

?x1?x2?4x3?x4?1??x1?x2?x3?4x4?1其精确解为x=[-1,-1,-1,-1]’.

解：令A=[-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]=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表示迭代次数上限，若迭代次数大于N，则迭代失败 %emg表示控制精度 %w表示松弛因子

%用SOR迭代法求线性方程组A*x=b的解 %k表示迭代次数

%x表示用迭代法求得的线性方程组的近似解 n=length(A);

18

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 j>=i

sum=sum+A(i,j)*x1(j); elseif j

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; end end x=x1;

算例15：用超松弛(SOR)迭代法求解方程组：

??4x1?x2?x3?x4?1?x?4x?x?x?1?1234 ?

?x1?x2?4x3?x4?1??x1?x2?x3?4x4?1其精确解为x=[-1,-1,-1,-1]’.

解：令A=[-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)

当松弛因子为1时，超松弛迭代法等同于Gauss-Seidel迭代法.

5.4 对称超松弛迭代

用对称超松弛(SSOR)迭代法求解线性方程组. MATLAB文件：（文件名：SSORmethod.m） function [x,k]=SSORmethod(A,b,x0,N,emg,w) %A是线性方程组的左端矩阵

19

%b是右端向量 %x0是迭代初始值

%N表示迭代次数上限，若迭代次数大于N，则迭代失败 %emg表示控制精度 %w表示松弛因子

%用SSOR迭代法求线性方程组A*x=b的解 %k表示迭代次数

%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

sum=sum+A(i,j)*x2(j); end end

x2(i)=(1-w)*x1(i)+w*(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

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; end end

20

b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0; end end

%回代求解

x(n)=b(n)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n

sum=sum+A(i,j)*x(j); end

x(i)=(b(i)-sum)/A(i,i); end

算例20：分别用Gauss列主元消去法和Gauss消去法求解线性方程组：

?0.3?10?15?5.291 ??11.2?1??1??x1??59.17???????6.13?12??x2??46.78??

952??x3??1??????x2211????4??59.143解：令A=[0.3*10^-15,59.14,3,1;5.291,-6.13,-1,2;11.2,9,5,2;1,2,1,1]; b=[59.17,46.78,1,2]’;

在命令窗口输入：

Gauss消去法: x=lu_decompose(A,b)

Gauss列主元消去法: x=Gauss_pivot(A,b)

二者的结果比较，用Gauss列主元消去法所得的结果要更准确.

26