matlab用于计算方法的源程序(5)

%x0:迭代初始值，列向量

%eps:误差限，可缺省，缺省值为0.5e-6 %

%应用举例:

%A=[10 3 1;2 -10 3;1 3 10];b=[14;-5;14];x0=[0;0;0]; %[x k]=EqtsJacobi(A,b,x0,0.5e-6) %x=EqtsJacobi(A,b,x0)

if nargin==3 eps=\end

%检查输入参数 n=length(b);

if size(A,1) ~= n || n ~= length(x0) disp('输入参数有误！'); x=' '; k=' '; return; end

%迭代求解 k=0;

x=zeros(n,1); while 1 k=\ for i=\ z=0; for j=\

z=\ end

for j=\

z=\ end

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

if norm(x-x0)<=eps || k>30 break; end x0=x; end

if k>=30

disp('迭代次数太多！')

x=' '; return; end return;

17、Guass-Seidel迭代法求解线性方程组

function [x k]=EqtsGS(A,b,x0,eps) %Guass-Seidel迭代法求解线性方程组Ax=b %[x k]=EqtsGS(A,b,x0,eps) %x:解向量，列向量 %k:迭代次数 %A:系数矩阵 %b:列向量

%x0:迭代初始值，列向量

%eps:误差限，可缺省，缺省值为0.5e-6 %

%应用举例:

%A=[10 3 1;2 -10 3;1 3 10];b=[14;-5;14];x0=[0;0;0]; %[x k]=EqtsGS(A,b,x0,0.5e-6) %x=EqtsGS(A,b,x0)

if nargin==3 eps=\end

%检查输入参数 n=length(b);

if size(A,1) ~= n || n ~= length(x0) disp('输入参数有误！'); x=' '; k=' '; return; end

%迭代求解 k=0;

x=zeros(n,1); while 1 k=\

for i=\ z=0; for j=\

z=\ end

for j=\

z=\ end

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

if norm(x-x0)<=eps || k>30 break; end x0=x; end

if k>=30

disp('迭代次数太多！') x=' '; return; end return;

ps:其实这段程序与Jacobi迭代法的程序相比，只修改了一个字母，不过，收敛速度却有明显提高

16、超松弛（SOR,Successive Over-Relaxation）迭代法求解线性方程组

function [x k]=EqtsSOR(A,b,x0,omiga,eps)

%超松弛（SOR,Successive Over-Relaxation）迭代法求解线性方程组Ax=b %[x k]=EqtsSOR(A,b,x0,eps) %x:解向量，列向量 %k:迭代次数 %A:系数矩阵 %b:列向量

%x0:迭代初始值，列向量

%omiga:松弛因子，可缺省，缺省值为1，即为GS迭代法 %eps:误差限，可缺省，缺省值为0.5e-6 %

%应用举例:

%A=[4 3 0;3 4 -1;0 -1 4];b=[24;30;-24];x0=[1;1;1];omiga=1.25; %[x k]=EqtsSOR(A,b,x0,omiga,0.5e-6) %x=EqtsSOR(A,b,x0)

if nargin==4 eps=\end

if nargin==3 omiga=\ eps=\end

%检查输入参数 n=length(b);

if size(A,1) ~= n || n ~= length(x0) disp('输入参数有误！'); x=' '; k=' '; return; end

%迭代求解 k=0;

x=zeros(n,1); while 1 k=\ for i=\ z=0; for j=\

z=\ end

for j=\

z=\ end

x(i)=(1-omiga)*x0(i)+omiga*(b(i)-z)/A(i,i); end

if norm(x-x0)<=eps || k>30 break; end x0=x; end

if k>=30

disp('迭代次数太多！')

x=' '; return; end return;

17、复化梯形公式求积分 function y=\%求积公式的测试函数，被积函数 y=3^x; return;

function y=\%求积公式被积函数的二次导数的相反值 y=-log(3)*log(3)*3^x; return;

function y=\%求积公式被积函数的四次导数的相反值 y=-(log(3))^4*3^x; return;

以上几个函数用于下列函数的测试

function [T n]=IntCompTrape(f,D2f,a,b,eps) %复化梯形公式求积分

%[T n]=IntCompTrape(f,D2f,a,b,eps) %T:求积结果 %n:区间等分数

%f:被积函数，可利用函数脚本文件事先定义，也可以利用内联函数

?f:被积函数的二次导数的相反值，可利用函数脚本文件事先定义，也可以利用内联函数 % 取相反值是为了便于计算被积函数的二次导数在区间[a,b]上的最大值 %a:积分下限 %b:积分上限