几种常用数值积分方法的比较(8)

贵州师范学院毕业论文（设计）

4 龙贝格积分方法的MATLAB实现 构造数表来逼近积分?af(x)dx?R(J,J)

其中。R(J,J)表示数表的最后一行，最后一列的值。 程序四：

function [R,quad,err,h]=romber(f,a,b,n,delta) % f是被积函数

% a,b分别是积分的上下限 % n+1是T数表的列数 % delta是允许误差 % R是T数表 % quad是所求积分值 M=1; h=b-a; err=1 J=0; R=zeros(4,4);

R(1,1)=h*(feval('f',a)+feval('f',b))/2 while ((err>delta)&(J

b24

贵州师范学院毕业论文（设计）

for p=1:M

x=a+h*(2*p-1); s=s+feval('f',x); end

R(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J

R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1); end

err=abs(R(J,J)-R(J+1,K+1));

end

quad=R(J+1,J+1)

先用M文件定义一个名为f.m的函数：

function y=f(x) if x==0 y=1;

else

y=sin(x)/x;

end

在MATLAB命令窗口中输入

>> romber('f',0,1,5,0.5*(10^(-8))) 回车得到 如图4.1

25

贵州师范学院毕业论文（设计）

图4.1

5 高斯-勒让德求积公式的MATLAB实现 程序五：

function [A,x]=Guass1(N) i=N+1;

f=((sym('t'))^2-1)^i; f=diff(f,i); t=solve(f); for j=1:i for k=1:i

X(j,k)=t(k)^(j-1); end

if mod(j,2)==0 B(j)=0;

26

贵州师范学院毕业论文（设计）

else B(j)=2/j; end end X=inv(X); for j=1:i A(j)=0; x(j)=0; for k=1:i

A(j)=A(j)+X(j,k)*B(k); x(j)=x(j)+t(j); end x(j)=x(j)/k; end

function g= GuassLegendre (a,b,n,m)

% a,b分别是积分的上下限； % n+1为节点个数；

% m是调用f1.m中第几个被积函数； [A,x]=Guass1(n); g=0; for i=1:n+1

y(i)=(b-a)/2*x(i)+(a+b)/2; f(i)=f1(m,y(i));

27

贵州师范学院毕业论文（设计）

g=g+(b-a)/2*f(i)*A(i); end

用M文件分别把上面两个自定义函数定义为名为Guass1.m函数和GuassLegendre.m函数

用M文件定义一个名为f1.m的函数function f=f1(i,x) g(1)=sqrt(x); if x==0 g(2)=1; else

g(2)=sin(x)/x; end

g(3)=4/(1+x^2); f=g(i);

在MATLAB命令窗口中输入

>> GuassLegendre (0,1,2,2) >> GuassLegendre (0,1,3,2) 回车得到 如图5.1

28

贵州师范学院毕业论文（设计）

图5.1

29

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库几种常用数值积分方法的比较(8)在线全文阅读。