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=\end tic; k=0; while 1

x=\ k=\

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

if k >= 30

disp('迭代次数太多。'); x=\ t=\end

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

function y=\

%牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量

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=\

%牛顿迭代法解非线性方程组的测试函数的导数 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=\end tic; k=0; while 1

x=\此处可采用其他方法避免求逆 k=\

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

if k >= 15

disp('迭代次数太多。'); x=\ t=\end

3、Lagrange插值法

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

function f=\%构造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=\else

disp('节点个数和函数值个数不同！') f=' '; return; end p=0; for i=\ l=\ for j=\ if j==i continue; end

%利用卷积计算Lagrange基函数 l=conv(l,[1 -x(j)]./(x(i)-x(j))); end

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

if nargin==3

f=\计算插值多项式在x0处的值 else

f=\把插值多项式的向量形式转化为插值多项式的符号形式 end

function f=\%构造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=\else

disp('节点个数和函数值个数不同！') f=' '; return; end syms t; f=0; for i=\ l=\ for j=\ if j==i continue; end

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

f=\

simplify(f);%化简多项式 if i==n

if nargin==3

f=\计算插值多项式f在点x0处的值 else

f=\计算插值多项式，展开并合并同类项 f=\设置多项式系数的有效数字 end end end

4、Newdon插值法

function f=\%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=\else

disp('节点个数和函数值个数不同！') f=' '; return; end

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

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

%计算差商矩阵

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

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=\

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

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

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

if nargin==3

f=\计算插值多项式在x0处的值 else

f=\把插值多项式的向量形式转化为插值多项式的符号形式 end

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

function x=\%基本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=\ for i=\

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

for i=\ for j=\

A(i,j)=A(i,j)-l(i,k)*A(k,j); end for j=\ A(i,j)=0; end end end

%回代过程 x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n); for i=\ y=\ for j=\ y=y+A(i,j)*x(j); end

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

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

function x=\%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=\ for j=\ z=0;

for r=\

z=\ end

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

for i=\ z=0;

for r=\

z=\ end

A(i,k)=(A(i,k)-z)/A(k,k); end end %求解

x=zeros(size(b)); for i=\ z=\ for k=\ z=z+A(i,k)*x(k); end x(i)=b(i)-z; end

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

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

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

function x=\%追赶法求解三对角线性方程组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=\

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=\

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

%赶的过程 x(n)=x(n)/D(n); for i=\

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

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

function x=\%主元素的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=\ %换行

[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=\ c=A(idx1,j); A(idx1,j)=A(idx2,j); A(idx2,j)=c; end

for i=\

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

for i=\ for j=\

A(i,j)=A(i,j)-l(i,k)*A(k,j); end for j=\ A(i,j)=0; end end end

%回代过程 x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n); for i=\ y=\ for j=\ y=y+A(i,j)*x(j); end

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

9、Euler法求解常微分方程

function outXY=\%简单欧拉法求解常微分方程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=\end

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

outXY(1,1)=x0; outXY(1,2)=y0; for i=\

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=\end tic;

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

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

k=\记录二分次数

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=\ else a=\ end end t=\ 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=\end tic; k=0; while 1

x=\ k=\

if abs(x-x1) < eps || k > 30 break; end x0=x1; x1=x; end t=toc;

if k >= 30

disp('迭代次数太多。'); x=\ t=\end

12、阻尼Newdon法求解非线性方程组 function y=\%带阻尼因子的牛顿迭代法测试方程组 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=\

%带阻尼因子的牛顿迭代法测试，方程组的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=\ eps=\end

if nargin==4 eps=\end

I=eye(size(df(x0))); tic; k=0; while 1

x=\此处可采用其他方法避免求逆 k=\

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

x0=x; end t=toc;

if k >= 30

disp('迭代次数太多。'); x=\ t=\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=\ eps=\end

if nargin==4 eps=\

end tic; k=0; while 1

x=\此处可采用其他方法避免求逆 k=\

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

if k >= 30

disp('迭代次数太多。'); x=\ t=\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=\end

x_const=x0; tic; k=0;

A=inv(df(x_const)); while 1

x=\此处可采用其他方法避免求逆 k=\

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

if k >= 30

disp('迭代次数太多。'); x=\ t=\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=\end tic; k=0;

B0=inv(df(x0)); while 1

x=\ k=\

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

B=\ x0=x; B0=B; end t=toc;

if k >= 30

disp('迭代次数太多。'); x=\ t=\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=\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:积分上限

%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=\默认精度 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=\ x1=a+k*h; T=\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=\默认精度 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=\ x1=a+k*h; x2=x1+h/2;

S=\end S=S*h/6; return;

本文来源：https://www.bwwdw.com/article/xab7.html