数值计算方法中的一些常用算法的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;