几种常用数值积分方法的比较 下载本文

贵州师范学院毕业论文(设计)

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