数值计算方法中的一些常用算法的Matlab源码，这些程序都是原创，传上来仅供大家参考，不足之处请大家指正，切勿做其它用途……

说明：这些程序都是脚本函数，不可直接运行，需要创建函数m文件，保存时文件名必须与函数名相同，懂一点儿Matlab的朋友应该知道。每个程序的说明里面都附了测试例子

1、Newdon迭代法求解非线性方程

function [x k t]=NewdonToEquation(f,df,x0,eps)

%牛顿迭代法解线性方程

%[x k t]=NewdonToEquation(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:原函数，定义为内联函数

�:函数的倒数，定义为内联函数

%x0:初始值

%eps:误差限

%

%应用举例：

%f=inline('x^3+4*x^2-10');

�=inline('3*x^2+8*x');

%x=NewdonToEquation(f,df,1,0.5e-6)

%[x k]=NewdonToEquation(f,df,1,0.5e-6)

%[x k t]=NewdonToEquation(f,df,1,0.5e-6)

%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6

%[x k t]=NewdonToEquation(f,df,1)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-f"(x0)./df(x0);

k="k"+1;

if abs(x-x0) <

eps || k >30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="0";

t="0";

end

2、Newdon迭代法求解非线性方程组

function y="NewdonF"(x)

%牛顿迭代法解非线性方程组的测试函数

%定义是必须定义为列向量

y(1,1)=x(1).^2-10*x(1)+x(2).^2+8;

y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8;

return;

function y="NewdonDF"(x)

%牛顿迭代法解非线性方程组的测试函数的导数

y(1,1)=2*x(1)-10;

y(1,2)=2*x(2);

y(2,1)=x(2).^+1;

y(2,2)=2*x(1).*x(2)-10;

return;

以上两个函数仅供下面程序的测试

function [x k t]=NewdonToEquations(f,df,x0,eps)

%牛顿迭代法解非线性方程组

%[x k t]=NewdonToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组(事先定义)

�:方程组的导数(事先定义)

%x0:初始值

%eps:误差限

%

%说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义

%  另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示

%

%应用举例：

%x0=[0,0];eps=0.5e-6;

%x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6

%[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-inv"(df(x0))*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0)

< eps || k > 15

break;

end

x0=x;

end

t=toc;

if k >= 15

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

3、Lagrange插值法

提供两个程序，采用了不同的方法

function f="InterpLagrange"(x,y,x0)

%构造Lagrange插值多项式

%此函数中借助向量卷积来求Lagrange基函数，运算速度较快

%f=InterpLagrange(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式:

%f=InterpLagrange(x,y) 返回插值多项式

%f=InterpLagrange(x,y,x0) 返回插值多项式在点x0处的值

%举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33;

%f=InterpLagrange(x,y)

%f=InterpLagrange(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

p=0;

for i="1:n"

l="y"(i);

for j="1:n"

if

j==i

continue;

end

%利用卷积计算Lagrange基函数

l=conv(l,[1

-x(j)]./(x(i)-x(j)));

end

%p是一向量，表示插值多项式的系数

p="p"+l;

end

if nargin==3

f="polyval"(p,x0);%计算插值多项式在x0处的值

else

f="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式

end

function f="InterpLagrange2"(x,y,x0)

%构造Lagrange插值多项式

%此函数中借助符号运算来求Lagrange基函数，运算速度较慢，不推荐此种方法

%f=InterpLagrange2(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式:

%f=InterpLagrange2(x,y) 返回插值多项式

%f=InterpLagrange2(x,y,x0) 返回插值多项式在点x0处的值

%举例:

%x=[0.32 0.34 0.36];y=[0.314567 0.333487 0.352274];x0=0.33;

%f=InterpLagrange2(x,y)

%f=InterpLagrange2(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

syms t;

f=0;

for i="1:n"

l="y"(i);

for j="1:n"

if

j==i

continue;

end

l=l*(t-x(j))/(x(i)-x(j));%借助符号运算，计算Lagrange基函数

end

f="f"+l;

simplify(f);%化简多项式

if i==n

if

nargin==3

f="subs"(f,'t',x0);%计算插值多项式f在点x0处的值

else

f="collect"(f);%计算插值多项式，展开并合并同类项

f="vpa"(f,6);%设置多项式系数的有效数字

end

end

end

4、Newdon插值法

function f="InterpNewdon"(x,y,x0)

%Newdon插值多项式

%f=InterpNewdon(x,y,x0)

%f:插值多项式或者是插值多项式在x0处的值

%x:节点

%y:函数值

%x0:某一测试点

%

%调用格式

%f=InterpNewdon(x,y) 返回插值多项式

%f=InterpNewdon(x,y,x0) 返回插值多项式在x0点的值

%应用举例:

%x=[1 2 3 4 5];y=[1 4 7 8 6];x0=6;

%f=InterpNewdon(x,y)

%f=InterpNewdon(x,y,x0)

if length(x)==length(y)

n="length"(x);

else

disp('节点个数和函数值个数不同！')

f=' ';

return;

end

A=zeros(n);%初始化差商矩阵

for i="1:n"

A(i,1)=y(i);%差商矩阵的第一列是函数值

end

%计算差商矩阵

%差商矩阵中对角线上的元素为Newdon插值多项式的系数

for j="2:n"

for i="j:n"

A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x(i)-x(i-j+1));

end

end

%求Newdon插值多项式

p=zeros(1,n);

for i="1:n"

p1=A(i,i);%差商矩阵对角线上的元素就是Newdon插值多项式的系数

for j="1:i-1"

p1=conv(p1,[1 -x(j)]);%计算Newdon插值多项式的基项

end

p1=[zeros(1,n-i),p1];%向量相加，维数必须相同。把向量的元素补齐

p="p"+p1;

end

if nargin==3

f="polyval"(p,x0);%计算插值多项式在x0处的值

else

f="poly2str"(p,'x');%把插值多项式的向量形式转化为插值多项式的符号形式

end

5、基本Guass消去法求解线性方程组

function x="EqtsBasicGuass"(A,b)

%基本Guass消去法求解线性方程组Ax=b

%x=EqtsBasicGuass(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[2 2 3;4 7 7;-2 4 5]; b=[3;1;-7];

%x=EqtsBasicGuass(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%(A|b)

A=[A b];

%消去过程

n=size(A,1);

l=zeros(n);

for k="1:n-1"

for i="k"+1:n

l(i,k)=A(i,k)/A(k,k);

end

for i="k"+1:n

for

j="k"+1:n+1

A(i,j)=A(i,j)-l(i,k)*A(k,j);

end

for

j="1:k"

A(i,j)=0;

end

end

end

%回代过程

x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n);

for i="n-1:-1:1"

y="0";

for j="i"+1:n

y=y+A(i,j)*x(j);

end

x(i)=(A(i,n+1)-y)/A(i,i);

end

return;

6、三角分解法求解线性方程组

function x="EqtsDoolittleLU"(A,b)

%Doolittle分解法求解线性方程组Ax=b

%x=EqtsDoolittleLU(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[6 2 1 -1;2 4 1 0;1 1 4 -1;-1 0 -1 3];b=[6;-1;5;-5];

%x=EqtsDoolittleLU(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%分解

%把L和U的元素存储在A的相应位置上

n=length(b);

for k="1:n"

for j="k:n"

z=0;

for

r="1:k-1"

z="z"+A(k,r)*A(r,j);

end

A(k,j)=A(k,j)-z;

end

for i="k"+1:n

z=0;

for

r="1:k-1"

z="z"+A(i,r)*A(r,k);

end

A(i,k)=(A(i,k)-z)/A(k,k);

end

end

%求解

x=zeros(size(b));

for i="1:n"

z="0";

for k="1:i-1"

z=z+A(i,k)*x(k);

end

x(i)=b(i)-z;

end

for i="n:-1:1"

z="0";

for k="i"+1:n

z=z+A(i,k)*x(k);

end

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

end

return

7、追赶法求解三对角线性方程组

function x="EqtsForwardAndBackward"(L,D,U,b)

%追赶法求解三对角线性方程组Ax=b

%x=EqtsForwardAndBackward(L,D,U,b)

%x:三对角线性方程组的解

%L:三对角矩阵的下对角线，行向量

%D:三对角矩阵的对角线，行向量

%U:三对角矩阵的上对角线，行向量

%b:线性方程组Ax=b中的b，列向量

%

%应用举例:

%L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]';

%x=EqtsForwardAndBackward(L,D,U,b)

%检查参数的输入是否正确

n=length(D);m=length(b);

n1=length(L);n2=length(U);

if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m

disp('输入参数有误！')

x=' ';

return;

end

%追的过程

for i="2:n"

L(i-1)=L(i-1)/D(i-1);

D(i)=D(i)-L(i-1)*U(i-1);

end

x=zeros(n,1);

x(1)=b(1);

for i="2:n"

x(i)=b(i)-L(i-1)*x(i-1);

end

%赶的过程

x(n)=x(n)/D(n);

for i="n-1:-1:1"

x(i)=(x(i)-U(i)*x(i+1))/D(i);

end

return;

8、主元素的Guass消去法求解线性方程组

function x="EqtsClmnPrimElemGuass"(A,b)

%主元素的Guass消去法求解线性方程组Ax=b

%x=EqtsClmnPrimElemGuass(A,b)

%x:解向量，列向量

%A:线性方程组的矩阵

%b:列向量

%

%应用举例:

%A=[-0.002 2 2;1 0.78125 0;3.996 5.5625 4];

%b=[0.4;1.3816;7.4178];

%x=EqtsClmnPrimElemGuass(A,b)

%检查输入参数

if size(A,1) ~= size(b,1)

disp('输入参数有误！');

x=' ';

return;

end

%(A|b)

A=[A b];

%消去过程

n=size(A,1);

l=zeros(n);

for k="1:n-1"

%换行

[a

idx1]=max(abs(A(k:n,k)));%寻找绝对值最大的元素的下标

[b

idx2]=min(abs(A(k:n,k)));%寻找绝对值最小的元素的下标

idx1=idx1+k-1;

idx2=idx2+k-1;

for j="1:n"+1

c=A(idx1,j);

A(idx1,j)=A(idx2,j);

A(idx2,j)=c;

end

for i="k"+1:n

l(i,k)=A(i,k)/A(k,k);

end

for i="k"+1:n

for

j="k"+1:n+1

A(i,j)=A(i,j)-l(i,k)*A(k,j);

end

for

j="1:k"

A(i,j)=0;

end

end

end

%回代过程

x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n);

for i="n-1:-1:1"

y="0";

for j="i"+1:n

y=y+A(i,j)*x(j);

end

x(i)=(A(i,n+1)-y)/A(i,i);

end

return;

9、Euler法求解常微分方程

function outXY="ODEEuler"(f,x0,y0,h,PointNum)

%简单欧拉法求解常微分方程dy/dx=f

%outXY=ODEEuler(f,x0,y0,h,PointNum)

%outXY:所取点的横纵坐标。第一列为横坐标，第二列为纵坐标

%f:函数f(x,y)，可利用脚本函数文件事先定义，也可利用内联函数

%x0:初始值的横坐标

%y0:初始值的纵坐标

%h:步长

%PointNum:计算步数，默认为30

%

%应用举例:

%f=inline('y-2*x/y','x','y');

%out=ODEEuler(f,0,1,0.1,10)

if nargin==4

PointNum="30";

end

outXY=zeros(PointNum+1,2);%初始化

outXY(1,1)=x0;

outXY(1,2)=y0;

for i="1:PointNum"

outXY(i+1,2)=outXY(i,2)+h*f(outXY(i,1),outXY(i,2));%简单Euler公式

outXY(i+1,1)=outXY(i,1)+h;

end

10、二分法解非线性方程

function [x k t]=DichotomyToEquation(f,a,b,eps)

%使用二分法解非线性方程

%[x k t]=DichotomyToEquation(f,a,b,eps)

%x:近似解

%k:二分次数

%t:运算时间

%f:函数，定义为内联函数

%a,b:区间端点

%eps:误差限

%

%应用举例:

%f=inline('x^3+4*x^2-10');

%x=DichotomyToEquation(f,1,2,0.5e-6)

%[x k]=DichotomyToEquation(f,1,2,0.5e-6)

%[x k t]=DichotomyToEquation(f,1,2,0.5e-6)

%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6

%[x k t]=DichotomyToEquation(f,1,2)

if nargin==3

eps="0".5e-6;

end

tic;

if f(a)*f(b) > 0

disp('区间太大或在此区间内无零点。');

else

k="0";%记录二分次数

while 1

x=(b+a)/2;

k=k+1;

if abs(x-a)

< eps || k > 30

break;

end

if f(a)*f(x)

< 0

b="x";

else

a="x";

end

end

t="toc";

x=(b+a)/2;

if k >=

30

disp('迭代次数太多。');

x=0;

t=0;

end

end

11、弦割法求解非线性方程

function [x k t]=ChordsecantToEquation(f,x0,x1,eps)

%弦割法求解非线性方程

%[x k t]=ChordsecantToEquation(f,x0,x1,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:原函数，定义为内联函数

%x0,x1:初始值

%eps:误差限

%

%应用举例：

%f=inline('x^3+4*x^2-10');

%x=ChordsecantToEquation(f,1,2,0.5e-6)

%[x k]=ChordsecantToEquation(f,1,2,0.5e-6)

%[x k t]=ChordsecantToEquation(f,1,2,0.5e-6)

%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6

%[x k t]=ChordsecantToEquation(f,1,2)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x1-f"(x1)*(x1-x0)./(f(x1)-f(x0));

k="k"+1;

if abs(x-x1) <

eps || k > 30

break;

end

x0=x1;

x1=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="0";

t="0";

end

12、阻尼Newdon法求解非线性方程组

function y="NewdonDampingF"(x)

%带阻尼因子的牛顿迭代法测试方程组

y(1,1)=x(1)^2-10*x(1)+x(2)^2+23;

y(2,1)=x(1)*x(2)^2+x(1)-10*x(2)+2;

return;

function y="NewdonDampingDF"(x)

%带阻尼因子的牛顿迭代法测试，方程组的Jacobi矩阵

y(1,1)=2*x(1)-10;

y(1,2)=2*x(2);

y(2,1)=x(2)^2+1;

y(2,2)=2*x(1)*x(2)-10;

return;

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

function [x k t]=NewdonDampingToEquations(f,df,x0,yita,eps)

%带阻尼因子的Newdon迭代格式求解非线性方程组

%[x k t]=NewdonDampingToEquations(f,df,x0,yita,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

�:方程组的Jacobi矩阵

%x0:初始值

%yita:阻尼因子，为了克服Jacobi矩阵的奇异性

%eps:误差限

%

%应用举例：

%x0=[2.5;2.5];yita=1e-5;eps=0.5e-6;

%x=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps)

%[x

k]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps)

%[x k

t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita,eps)

%函数的最后两个参数也可以不写，默认情况下，yita=1e-4;eps=0.5e-6

%[x k

t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0)

%[x k

t]=NewdonDampingToEquations(@NewdonDampingF,@NewdonDampingDF,x0,yita)

if nargin==3

yita="1e-4";

eps="0".5e-6;

end

if nargin==4

eps="0".5e-6;

end

I=eye(size(df(x0)));

tic;

k=0;

while 1

x="x0-inv"(df(x0)+yita*I)*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0)

< eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

13、牛顿下山法求解非线性方程组

测试函数见2(Newdon迭代法求解非线性方程组)

function [x k t]=NewdonDescendToEquations(f,df,x0,omiga,eps)

%下降牛顿迭代法求解非线性方程组

%[x k t]=NewdonDescendToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

�:方程组的Jacobi矩阵

%x0:初始值

%omiga:下降因子，通常在区间(0,1)内选择。当omiga=1时，即为Newdon迭代格式

%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps)

%[x

k]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps)

%[x k

t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga,eps)

%函数的最后两个参数也可以不写，默认情况下，omiga=0.8;eps=0.5e-6

%[x k t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0)

%[x k

t]=NewdonDescendToEquations(@NewdonF,@NewdonDF,x0,omiga)

if nargin==3

omiga="0".8;

eps="0".5e-6;

end

if nargin==4

eps="0".5e-6;

end

tic;

k=0;

while 1

x="x0-omiga"*inv(df(x0))*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0)

< eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

14、简化牛顿法求解非线性方程组

测试函数见2(Newdon迭代法求解非线性方程组)

function [x k t]=NewdonSimplifyToEquations(f,df,x0,eps)

%简化牛顿格式求解非线性方程组

%[x k t]=NewdonSimplifyToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

�:方程组的Jacobi矩阵

%x0:初始值

%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k t]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0,eps)

%函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6

%[x k t]=NewdonSimplifyToEquations(@NewdonF,@NewdonDF,x0)

if nargin==3

eps="0".5e-6;

end

x_const=x0;

tic;

k=0;

A=inv(df(x_const));

while 1

x="x0-A"*f(x0);%此处可采用其他方法避免求逆

k="k"+1;

if norm(x-x0)

< eps || k > 30

break;

end

x0=x;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

15、逆Broyden秩1方法求解非线性方程组

测试函数见2(Newdon迭代法求解非线性方程组)

function [x k t]=Broyden1InvToEquations(f,df,x0,eps)

%逆Broyden秩1方法求解非线性方程组

%function [x k t]=Broyden1InvToEquations(f,df,x0,eps)

%x:近似解

%k:迭代次数

%t:运算时间

%f:方程组

�:方程组的Jacobi矩阵

%x0:初始值

%eps:误差限

%

%应用举例：

%x0=[0;0];eps=0.5e-6;

%x=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps)

%[x k t]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0,eps)

%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6

%[x k t]=Broyden1InvToEquations(@NewdonF,@NewdonDF,x0)

if nargin==3

eps="0".5e-6;

end

tic;

k=0;

B0=inv(df(x0));

while 1

x="x0-B0"*f(x0);

k="k"+1;

if norm(x-x0)

< eps || k> 30

break;

end

s="x-x0";

y="f"(x)-f(x0);

B="B0"+(s-B0*y)*s'*B0/(s'*B0*y);

x0=x;

B0=B;

end

t=toc;

if k >= 30

disp('迭代次数太多。');

x="zeros"(size(x0));

t="0";

end

16、Jacobi迭代法求解线性方程组

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

%Jacobi迭代法求解线性方程组Ax=b

%[x k]=EqtsJacobi(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]=EqtsJacobi(A,b,x0,0.5e-6)

%x=EqtsJacobi(A,b,x0)

if nargin==3

eps="0".5e-6;

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="k"+1;

for i="1:n"

z=0;

for

j="1:i-1"

z="z"+A(i,j)*x0(j);

end

for

j="i"+1:n

z="z"+A(i,j)*x0(j);

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="0".5e-6;

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="k"+1;

for i="1:n"

z=0;

for

j="1:i-1"

z="z"+A(i,j)*x(j);

end

for

j="i"+1:n

z="z"+A(i,j)*x0(j);

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="0".5e-6;

end

if nargin==3

omiga="1";

eps="0".5e-6;

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="k"+1;

for i="1:n"

z=0;

for

j="1:i-1"

z="z"+A(i,j)*x(j);

end

for

j="i"+1:n

z="z"+A(i,j)*x0(j);

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="IntF"(x)

%求积公式的测试函数，被积函数

y=3^x;

return;

function y="IntD2F"(x)

%求积公式被积函数的二次导数的相反值

y=-log(3)*log(3)*3^x;

return;

function y="IntD4F"(x)

%求积公式被积函数的四次导数的相反值

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:积分上限

%eps:误差限

%

%应用举例:

%事先定义

%[T n]=IntCompTrape(@IntF,@IntD2F,0,1)

%[T n]=IntCompTrape(@IntF,@IntD2F,0,1,0.5e-7)

%利用内联函数

%F=inline('3^x');D2F=inline('-(log(3))^2*3^x');

%[T n]=IntCompTrape(F,D2F,0,1)

if nargin==4

eps="0".5e-7;%默认精度

end

%求被积函数的二次导数在区间[a,b]上的最大值

[x,fval]=fminbnd(D2f,a,b,optimset('TolX',eps/10));

fmax=-fval;

%计算等分区间数

n=ceil(sqrt(abs(fmax*(b-a)^3/12/eps)));

h=(b-a)/n;%步长

T=f(a)+f(b);

for k="1:n-1"

x1=a+k*h;

T="T"+2*f(x1);

end

T=h*T/2;

return;

18、复化Simpson公式求积分

测试函数见17

function [S n]=IntCompSimpson(f,D4f,a,b,eps)

%复化辛普森公式求积分

%[S n]=IntCompSimpson(f,D4f,a,b,eps)

%S:数值求积结果

%n:区间等分数

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

�f:被积函数的四次导数的相反值，可利用函数脚本文件事先定义，也可以利用内联函数

%  取相反值是为了便于计算被积函数的四次导数在区间[a,b]上的最大值

%a:积分下限

%b:积分上限

%eps:误差限

%

%应用举例:

%事先定义

%[T n]=IntCompSimpson(@IntF,@IntD4F,0,1)

%[S n]=IntCompSimpson(@IntF,@IntD4F,0,1,0.5e-7)

%利用内联函数

%F=inline('3^x');D4F=inline('-(log(3))^4*3^x');

%[S n]=IntCompSimpson(F,D4F,0,1)

if nargin==4

eps="0".5e-7;%默认精度

end

%求被积函数的四次导数在区间[a,b]上的最大值

[x,fval]=fminbnd(D4f,a,b,optimset('TolX',eps/10));

fmax=-fval;

%计算等分区间数

n=ceil(sqrt(sqrt(abs((b-a)^5*fmax/16/180/eps))));

h=(b-a)/n;%步长

S=f(a)+f(b)+4*f(a+h/2);

for k="1:n-1"

x1=a+k*h;

x2=x1+h/2;

S="S"+4*f(x2)+2*f(x1);

end

S=S*h/6;

return;