http://www.mathworks.com/matlabcentral/fileexchange/17648这里有其matlab代码实现，楼主参考下！
1. 用雅克比迭代法和高斯--赛德尔迭代法求解下列方程组，取迭代初值[000]。(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代法的收敛性，若均收敛，对比两种方法的收敛速度。
解：源程序：
①雅克比迭代法：建立函数文件jacobi.m
function [n,x]=jacobi(A,b,X,nm,w)
%用雅克比迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D)*(L+U) %计算迭代矩阵
g=inv(D)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
②高斯赛德尔迭代法：建立函数文件gaussseidel.m
function [n,x]=gaussseidel(A,b,X,nm,w)
%用高斯-赛德尔迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
I=eye(m) %生成m*m阶的单位矩阵
D=diag(diag(A))%令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-L)*U %计算迭代矩阵
g=inv(I-inv(D)*L)*(inv(D)*b) %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g %用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
运行过程及结果：
①雅克比迭代法的运行过程及结果为：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
X=[000]
nm=50
w=10^-6
>>jacobi(A,b,X,nm,w)
迭代次数为
n =
14
方程组的解为
x =
0.0254
0.5144
0.2026
②高斯赛德尔迭代法的运行过程及结果为：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
X=[000]
nm=50
w=10^-6
>>gaussseidel(A,b,X,nm,w)
迭代次数为
n =
6
方程组的解为
x =
0.0254
0.5144
0.2026
③用matlab中的A\b命令的运行过程及结果如下：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
>>A\b
ans =
0.0254
0.5144
0.2026
(1)由运行结果可知，编程得到的方程组的解为[0.02540.51440.2026]。而用matlab中的A\b命令求出的方程组的解为[0.02540.51440.2026]，完全一致。
(2)由运行过程可知，两种方法均收敛，雅克比迭代次数为14，高斯赛德尔迭代次数为6，说明后者的迭代速度比前者快。
2.用超松弛迭代法求解方程组，分别取松弛因子 ，取迭代初值[000]，迭代30次或满足 时停止计算。
(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代是否收敛，若收敛，松弛因子取何值时收敛最快，与有关定理的结论对照，看结果是否一致。
解：源程序：
①超松弛迭代法：建立函数文件sor22.m
function [n,x]=sor22(A,b,X,nm,w,ww)
%用超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-ww*L)*((1-ww)*D+ww*U) %计算迭代矩阵
g=ww*inv(D-ww*L)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
②向后超松弛迭代法：建立函数文件sorxh22.m
function [n,x]=sorxh22(A,b,X,nm,w,ww)
%用向后超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-ww*U)*((1-ww)*D+ww*L) %计算矩阵迭代矩阵
g=ww*inv(D-ww*U)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
③对称超松弛迭代法：建立函数文件ssor22.m
function [n,x]=ssor22(A,b,X,nm,w,ww)
%用对称超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
I=eye(m) %生成m*m阶的单位矩阵
D=diag(diag(A))%令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
%下面计算迭代矩阵和迭代格式中的常数项
M=inv(D-ww*U)*((1-ww)*D+ww*L)*inv(D-ww*L)*((1-ww)*D+ww*U)g=ww*inv(D-ww*U)*(I+((1-ww)*D+ww*L)*inv(D-ww*L))*b
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
由向后欧拉公式有:y(k+1)=y(k)-30*h*y(k+1)
变形求得:
y(k+1)=y(k)/(30*h+1)
故MATLAB程序有:
h=0.05
x=[0:h:1]
y(1)=1
for
k=1:length(x)-1
y(k+1)=y(k)/(30*h+1)
end
plot(x,y,'r.-')
title('向后欧拉')
grid
on
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）