Jacobi迭代法 Gauss-Seidel迭代法

Matlab线性方程组的迭代解法（Jacobi迭代法Gauss-Seidel迭代法）实验报告2008年11月09日星期日12:49
1.熟悉Jacobi迭代法，并编写Matlab程序matlab程序
按照算法(Jacobi迭代法)编写Matlab程序(Jacobi.m)
function [x, k, index]=Jacobi(A, b, ep, it_max)
%求解线性方程组的Jacobi迭代法，其中
% A ---方程组的系数矩阵
% b ---方程组的右端项
% ep ---精度要求。

省缺为1e-5
% it_max ---最大迭代次数，省缺为100
% x ---方程组的解
% k ---迭代次数
% index --- index=1表示迭代收敛到指定要求；
% index=0表示迭代失败
if nargin <4 it_max=100; end
if nargin <3 ep=1e-5; end
n=length(A); k=0;
x=zeros(n,1); y=zeros(n,1); index=1;
while 1
for i=1:n
y(i)=b(i);
for j=1:n
if j~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
if abs(A(i,i))<1e-10 | k==it_max
index=0; return;
end
y(i)=y(i)/A(i,i);
end
if norm(y-x,inf)<ep
break;
end
x=y; k=k+1;
end
用Jacobi迭代法求方程组
的解。

输入：
A=[4 3 0;3 3 -1;0 -1 4]；
b=[24;30;-24]；
[x, k, index]=Jacobi(A, b, 1e-5, 100)
输出：
x =
-2.9998
11.9987
-3.0001
k =
100
index =
2.熟悉Gauss-Seidel迭代法，并编写Matlab程序
function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp)
%Gauss-Seidel迭代法求解线性方程组
%A-系数矩阵b-右端向量x0-初始迭代点errorBound-近似精度maxSp-最大迭代次数
%v-近似解sN-迭代次数vChain-迭代过程的所有值
step=0;
error=inf;
s=size(A);
D=zeros(s(1));
vChain=zeros(15,3);%最多能记录15次迭代次数
k=1;
fx0=x0;
for i=1:s(1)
D(i,i)=A(i,i);
end;
L=-tril(A,-1);
U=-triu(A,1);
while error>=errorBound & step<maxSp
x0=inv(D)*(L+U)*x0+inv(D)*b;
vChain(k,:)=x0';
k=k+1;
error=norm(x0-fx0);
fx0=x0;
step=step+1;
end
v=x0;
sN=step;
用Gauss-Seidel迭代法求解上题的线性方程组，取。

输入：
A=[4 3 0;3 3 -1;0 -1 4]；
b=[24;30;-24]；
x0=[0;0;0]；
[v,sN,vChain]=gaussSeidel(A,b,x0,0.00001,11)
输出：
v =
0.6169
11.1962
-4.2056
sN =
11
vChain =
6.0000 10.0000 -6.0000 -1.5000 2.0000 -3.5000
4.5000 10.3333 -
5.5000 -1.7500 3.6667 -3.4167
3.2500 10.6111 -5.0833 -1.9583 5.0556 -3.3472
2.2083 10.8426 -4.7361 -2.1319 6.2130 -
3.2894
1.3403 11.0355 -4.4468 -
2.2766 7.1775 -
3.2411
0.6169 11.1962 -4.2056
0 0 0
0 0 0
0 0 0
0 0 0
s
数值实验
数值实验要求：
数值实验报告内容：要包含题目、算法公式、完整的程序、正确的数值结果和图形以及相应的误差分析。

在本课程网站上提交数值实验报告的电子文档。

一、为了逼近飞行中的野鸭的顶部轮廓曲线，已经沿着这条曲线选择了一组点。

见下表。

1．对这些数据构造三次自然样条插值函数，并画出得到的三次自然样条插值曲线；2．对这些数据构造Lagrang插值多项式，并画出得到的Lagrang插值多项式曲线。

x 0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0
f(x) 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25
x 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3
f(x) 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25
1．使用三次样条插值函数csape（）求解。

解：输入：
>>x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3]; >> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25]; >> pp=csape(x,y,'variational',[0 0])
得到：
pp =
form: 'pp'
breaks: [0.9000 1.3000 1.9000 2.1000 2.6000 3 3.9000 4.4000 4.7000 5 6 7 8 9.2000 10.5000 11.3000 11.6000 12 12.6000 13 13.3000]
coefs: [20x4 double]
pieces: 20
order: 4
dim: 1
再输入：
>> pp.coefs
得到：
ans =
-0.2476 0 0.5396 1.3000
0.9469 -0.2972 0.4208 1.5000
-2.9564 1.4073 1.0868 1.8500
-0.4466 -0.3666 1.2949 2.1000
0.4451 -1.0365 0.5934 2.6000
0.1742 -0.5025 -0.0222 2.7000
0.0781 -0.0322 -0.5034 2.4000
1.3142 0.0849 -0.4771
2.1500
-1.5812 1.2676 -0.0713 2.0500
0.0431 -0.1555 0.2623 2.1000
-0.0047 -0.0261 0.0808 2.2500
-0.0244 -0.0401 0.0146 2.3000
0.0175 -0.1135 -0.1390 2.2500
-0.0127 -0.0506 -0.3358 1.9500
-0.0203 -0.1002 -0.5318 1.4000
1.2134 -0.1490 -0.7312 0.9000
-0.8393 0.9431 -0.4929 0.7000
0.0364 -0.0640 -0.1413 0.6000
-0.4480 0.0014 -0.1789 0.5000
0.5957 -0.5361 -0.3928 0.4000
因此，三次样条函数S（x）为
最后输入：
>> xx=0:0.01:14;yy=interp1(x,y,xx,'csape');
>> plot(xx,yy)；xlabel('x');ylabel('y');
得到图形：
grange插值算法：
1) 输入N个节点数组；
2) 定义初始值和；
3)利用公式
求N 次插值基函数；
4)利用多项式求插值函数;
解：输入：
>>x=[0.9,1.3,1.9,2.1,2.6,3.0,3.9,4.4,4.7,5.0,6.0,7.0,8.0,9.2,10.5,11.3,11.6,12.0,12.6,13.0,13.3];
y=[1.3,1.5,1.85,2.1,2.6,2.7,2.4,2.15,2.05,2.1,2.25,2.3,2.25,1.95,1.4,0.9,0.7,0.6,0.5,0.4,0.25];
>>a=polyfit(x,y,length(x)-1);
>>p=vpa(poly2sym(a),5)
得到数值结果：p= .30732e-10*x^20+.42770e-8*x^19-.27708e-6*x^18+.11098e-4*x^17-.30784e-3*x^16+.62785 e-2*x^15-.97558e-1*x^14+1.1810*x^13-11.297*x^12+86.085*x^11-524.68*x^10+2558.0*x^9-9 942.3*x^8+30586.*x^7-73622.*x^6+.13627e6*x^5-.18907e6*x^4+.18913e6*x^3-.12803e6*x^2 +52170.*x-9593.4
再输入：
>> y1=polyval(a,x);
plot(x,y1);xlabel('x');ylabel('y');
得到图形：
结果分析：
由以上两图形可以看出三次样条插值的图形较lagrange插值的图形要光滑的多。

因为样条函数插值不仅具有较好的收敛性和稳定性，而且其光滑性也较高，因此样条函数成为了重要的插值工具，其中应用较多的是三次样条插值。

二、对于问题
将h=0.025的Euler法，h=0.05的改进的Euler法和h=0.1的4阶经典的Runge-Kutta法在这些方法的公共节点0.1，0.2，0.3，0.4和0.5处进行比较。

精确解为：。

1. Euler法
在x, y平面上微分方程的解在曲线上一点(x, y) 的切线斜率等于函数的值。

该曲线的顶点设为p ,再推进到p ( ),显然两个顶点p , p 的坐标有以下关系
Matlab程序：
function [x,y]=Euler(ydot,x0,y0,h,N)
% ydot为一阶微分方程的函数
% x0，y0为初始条件
% h为区间步长
% N为区间个数
% x为Xn构成的向量，y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(ydot,x(n),y(n));
end
解：取h=0.025，N=20，输入：
>> ydot=inline('y-x^2+1','x','y');
[t,u]=Euler(ydot,0,0.5,0.025,20)
得到数值结果：
t =
Columns 1 through 17
0 0.0250 0.0500 0.0750 0.1000 0.1250 0.1500 0.1750 0.2000 0.2250 0.2500 0.2750 0.3000 0.3250 0.3500 0.3750 0.4000
Columns 18 through 21
0.4250 0.4500 0.4750 0.5000
u =
Columns 1 through 17
0.5000 0.5375 0.5759 0.6153 0.6555 0.6966 0.7387 0.7816
0.8253 0.8700 0.9155 0.9618 1.0089 1.0569 1.1057 1.1553
1.2056
Columns 18 through 21
1.2568 1.3087 1.3613 1.4147
即采用Euler法得到：
u(0.1)=0.6555,u(0.2)=0.8253,u(0.3)=1.0089,u(0.4)=1.2056, u(0.5)=1.4147
2. 改进Euler方法
改进Euler公式.
y = yn+h/2(f ( ) + f (xn+1, + h* f ( ))) 即迭代公式为:
Matlab程序：
function [x,y]=Euler_ad(ydot,x0,y0,h,N)
% 改进Euler公式
% ydot为一阶微分方程的函数
% x0，y0为初始条件
% h为区间步长
% N为区间个数
% x为Xn构成的向量，y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
ybar=y(n)+h*feval(ydot,x(n),y(n));
y(n+1)=y(n)+h/2*(feval(ydot,x(n),y(n))+feval(ydot,x(n+1),ybar));
end
解：取h=0.05，N=10，输入：
>> ydot=inline('y-x^2+1','x','y');
[t,u]=Euler_ad(ydot,0,0.5,0.05,10)
得到数值结果：
t =
0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000
u =
0.5000 0.5768 0.6573 0.7414 0.8291 0.9202 1.0147 1.1126
1.2136 1.3178 1.4250
即采用改进的Euler法得到：
u(0.1)=0.6573,u(0.2)=0.8291,u(0.3)=1.0147,u(0.4)=1.2136,u(0.5)=1.4250
3.四阶Runge_Kutta法
求解公式为：
Matlab程序：
function [x,y]=Runge_Kutta4(ydot,x0,y0,h,N)
% 标准四阶Runge_Kutta方法
% ydot为一阶微分方程的函数
% x0，y0为初始条件
% h为区间步长
% N为区间个数
% x为Xn构成的向量，y为Yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(ydot,x(n),y(n));
k2=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k1);
k3=h*feval(ydot,x(n)+1/2*h,y(n)+1/2*k2);
k4=h*feval(ydot,x(n)+h,y(n)+k3);
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end
解：取h=0.1，N=5，输入：
>> ydot=inline('y-x^2+1','x','y');
[t,u]=Runge_Kutta4(ydot,0,0.5,0.1,5)
得到数值结果：
t =
0 0.1000 0.2000 0.3000 0.4000 0.5000
u =
0.5000 0.6574 0.8293 1.0151 1.2141 1.4256
结果比较：
t Euler法改进Euler法四阶Runge_Kutta 精确解
0.1 0.6555 0.6573 0.6574 0.6574
0.2 0.8253 0.8291 0.8293 0.8293
0.3 1.0089 1.0147 1.0151 1.0151
0.4 1.2056 1.2136 1.2141 1.2141
0.5 1.4147 1.4250 1.4256 1.4256
由以上结果可以看出改进的Euler法较Euler法计算精度有所提高，但还不是十分精确。

四
阶Runge_Kutta法具有非常高的精度，事实上，在求解微分方程初值问题，四阶Runge_Kutta 法是单步长中最优秀的方法，通常都用该方法求解实际问题。

三、
用Newton迭代法求方程的根时，分别取初始值，；
用Newton迭代法求方程时，分别取初始值，；
算法：
（1）取初始点x0 最大迭代次数N和精度要求ε，k=0.
（2）如果f’（xk）=0，则停止计算；否则计算
Xk+1=xk- f（xk）/f’（xk）
（3）若|xk+1- xk|<ε，则停止计算。

（4）若k=N，则停止计算；否则置k=k+1，转（2）。

Matlab程序：
function [x_star,index,it]= Newton(fun,x,ep,it_max)
% 求解非线性方程的Newton法
% fun（x) 为需要求根的函数，第一个分量是函数值，第二个分量是导数值
% fun=inline('[x^3-x-1,3*x^2-1]') 当f(x)=x^3-x-1;
% x为初始点
% ep为精度，缺省值为1e-5
% it_max为最大迭代次数，缺省值100
% x_star为当迭代成功时输出方程的根，失败时输出最后的迭代值
% index为指标变量，index=1表示迭代成功index=0表示失败
% it为迭代次数
if nargin<4 it_max=100;end
if nargin<3 ep=1e-5;end
index=0;k=1;
while k<=it_max
x1=x; f=feval(fun,x);
if abs(f(2))<ep break; end
x=x-f(1)/f(2);
if abs(x-x1)<ep index=1;break;end
k=k+1;
end
x_star=x;it=k;
解：（1）由于f(x)=arctan(x) , f’(x)= 1/1+x2 , 取初始值x0=1.45,输入
>> fun=inline('[atan(x),1/(1+x^2)]');
>> [x_star,index,it]=Newton(fun,1.45)
得到数值结果：
x_star =1.6281e+004
index = 0
it =7
取初始值x0=0.5,输入
>> fun=inline('[atan(x),1/(1+x^2)]');
>> [x_star,index,it]=Newton(fun,0.5)
得到数值结果：
x_star =0
index = 1
it =4
输入x=-3:0.05:3; y=atan(x);
plot(x,y);xlabel('x');ylabel('y');
得到图形：
(2) 由于f(x)=x3-x-3=0, f’(x)= 3x2-1 , 取初始值x0=0.0,输入
>> fun=inline('[x^3-x-3,3*x^2-1]');
>> [x_star,index,it]=Newton(fun,0.0)
得到数值结果：
x_star = -0.0074
index = 0
it =101
取初始值x0=2.0,输入
>> fun=inline('[x^3-x-3,3*x^2-1]');
>> [x_star,index,it]=Newton(fun,2.0)
得到数值结果：
x_star = 1.6717
index = 1
it =4
输入x=0:0.05:3; y=x.^3-x-3;
plot(x,y);xlabel('x');ylabel('y');
得到图形：
结果分析：
从以上结果可以看出初值的选取与Newton法的收敛很有关系。

初值选的不好，Nexton法将不收敛，它的收敛性是在跟a附近讨论的。

选取初值时一定要十分小心。

四、用Jacobi迭代和SOR法分别求解线性方程组（教科书第86页算例2），并验证、输出SOR法的松弛因子w和对应的迭代次数。

1. Jacobi 迭代法
Jacobi迭代法的算法为:
Matlab程序：
function[x,k,index]=Jacobi(A,b,ep,it_max)
% 求解线性方程组的Jacobi迭代法
% A为系数矩阵
% b为方程组右端项
% ep为精度要求，缺省值1e-5
% it_max为最大迭代次数，缺省值100
% x为方程组的解
% k为迭代次数
% index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。

if nargin<4 it_max=100;end
if nargin<3 ep=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);index=1;
while 1
for i=1:n
y(i)=b(i);
for j=1:n
if j~=i
y(i)=y(i)-A(i,j)*x(j);
end
end
if abs(A(i,i))<1e-10|k==it_max
index=0;return;
end
y(i)=y(i)/A(i,i);
end
if norm(y-x,inf)<ep
break;
end
x=y;k=k+1;
end
function [A,b]=shuru
% 自己定义输入矩阵A和向量b的函数
n=15;
for i=1:n
if i==1|i==15 z(i)=3;
else z(i)=2;
end
for j=1:n
if j==i A(i,j)=4;
elseif j==i+1|j==i-1 A(i,j)=-1;
else A(i,j)=0;
end
end
end
b=z';
解：输入：
>> [A,b]=shuru;ep=1e-6;
>> [x,k,index]=Jacobi(A,b,ep)
得到数值结果：
x =(1.0000,1.0000,…,1.0000)T
k = 19
index =1
2. sor法
SOR迭代法的算法为：
Matlab程序：
function [x,k,index]=SOR1(A,b,ep,w,it_max)
% 求解线性方程组的SOR迭代法
% A为系数矩阵
% b为方程组右端项
% ep为精度要求，缺省值1e-5
% w为超松弛因子，缺省值为1；
% it_max为最大迭代次数，缺省值100
% x为方程组的解
% k为迭代次数
% index为指标变量index=1表示迭代收敛到指定要求index=0表示迭代失败。

if nargin<5 it_max=100;end
if nargin<4 w=1;end
if nargin<3 ep=1e-5;end
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);index=1;
while 1
y=x;
for i=1:n
z=b(i);
for j=1:n
if j~=i
z=z-A(i,j)*x(j);
end
end
if abs(A(i,i))<1e-10|k==it_max
index=0;return;
end
z=z/A(i,i);x(i)=(1-w)*x(i)+w*z;
end
if norm(y-x,inf)<ep
break;
end
k=k+1;
end
解：取w=1.1，输入：
>> [A,b]=shuru;ep=1e-6;w=1.1;
>> [x,k,index]=sor(A,b,ep,w)
得到数值结果：
x =(1.0000,1.0000,…,1.0000)T
k = 11
index =1
依次取w=0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5 得到下表：
松弛因子w 迭代次数
0.8 19
0.9 16
1.0 13
1.1 11
1.2 14
1.3 17
1.4 23
1.5 27
结果分析：
从以上结果可以看出在求解相同问题时，sor法的迭代次数要比Jacobi迭代法少很多，计算量小许多。

此外可以看出松弛因子w的选取对迭代次数的影响也十分大。

在实际计算时，最优松弛因子很难事先确定，一般可用试算法取近似最优值。

function X=jacobi(A,b,P,delta,max1) % A是n维非奇异阵% B是n维向量% P是初值% delta是误差界% X为所求的方程组AX=B的近似解N=length(b); for k=1:max1 for j=1:N X(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j); end err=abs(norm(X'-P)); P=X'; if (err<delta) break end end X=X';k,err >> A=[4,1,-1;1,-5,-1;2,-1,-6]>> b=[13;-8;-2]>> P=[0;0;0]>> X=jacobi(A,b,P,10^(-4),20) k = 9 err = 2.5713e-005 X = 3.0000 2.0000 1.0000
nargin是用来判断输入变量个数的函数，这样就可以针对不同的情况执行不同的功能。

通常可以用他来设定一些默认值，如下面的函数。

例子，函数test1的功能是输出a和b的和。

如果只输入一个变量，则认为另一个变量为0，如果两个变量都没有输入，则默认两者均为0。

function y=test1(a,b)if nargin==0 a=0;b=0;elseif nargin==1 b=0;endy=a+b; s
（2）迭代解法
迭代解法非常适合求解大型系数矩阵的方程组。

在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。

i．Jacobi迭代法
对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：
x=D-1(L+U)x+D-1b
与之对应的迭代公式为：
x(k+1)=D-1(L+U)x(k)+D-1b
这就是Jacobi迭代公式。

如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。

Jacobi迭代法的MA TLAB函数文件Jacobi.m如下：
function [y,n]=jacobi(A,b,x0,eps)
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=B*x0+f;
n=n+1;
end
例用Jacobi迭代法求解下列线性方程组。

设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件Jacobi.m，命令如下：
A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)
ii．Gauss-Serdel迭代法
在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：
x(k+1)=(D-L)-1Ux(k)+(D-L)-1b
该式即为Gauss-Serdel迭代公式。

和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。

Gauss-Serdel迭代法的MA TLAB函数文件gauseidel.m如下：
function [y,n]=gauseidel(A,b,x0,eps)
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=G*x0+f;
n=n+1;
end
例用Gauss-Serdel迭代法求解下列线性方程组。

设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件gauseidel.m，命令如下：
A=[10,-1,0;-1,10,-2;0,-2,10];
b=[9,7,6]';
[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)
例分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。

命令如下：
a=[1,2,-2;1,1,1;2,2,1];
b=[9;7;6];
[x,n]=jacobi(a,b,[0;0;0])
[x,n]=gauseidel(a,b,[0;0;0])。