龙贝格算法求解椭球周长

数值积分：

在实际应用中经常应用到计算方法去求解一些不易测量的零件的周长或面积。已知一个椭圆形边框如下图所示，试用龙贝格算法求解这个边框的周长，要求结果精确到6 位有效数字。

3.1 数学原理：

龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。

（1）将区间[a,b] 划分为 n等分，分点： ；根据梯形公式，求出 ，再根据和之间的递推公式求出；

（2）设m 为加速次数， k为划分区间次数，则由加速公式： （k=1,2,... ）求出第 k次划分，第 m次加速次数的梯形值，这样不断地循环，直到求出在满足精度条件下的某个作为积分值为止。

题中要求的椭圆周长可用函数 来求得。

3.2 程序设计：

%主函数程序

a = 0;

b = 3.1415926/2;

epsilon = 0.5e-6;

f = @(x)sqrt(1-7/16.*sin(x).*sin(x));

y =32*Romberg(f,a,b,epsilon);

%后面是画出函数图像，注，不是积分函数图像。是被积函数图像

x = 0:0.01:3.1415926/2;

z = sqrt(1-7/16.*sin(x).*sin(x));

plot(x,z),xlabel('x'),ylabel('y'),title(['Result = ',num2str(y)])

 

%定义龙贝格函数子程序

function [R,k,T]=Romberg(fun,a,b,tol)

k=0; % 迭代次数

n=1; % 区间划分个数

h=b-a;

T=h/2*(fun(a)+fun(b));

err=1;

while err>=tol

    k=k+1;

    h=h/2;

    tmp=0;

    for i=1:n

        tmp=tmp+fun(a+(2*i-1)*h);

    end

    T(k+1,1)=T(k)/2+h*tmp;

    for j=1:k

        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);

        T1(k+1,j+1)=32*T(k+1,j+1)

    end

    n=n*2;

    err=abs(T(k+1,k+1)-T(k,k));

end

R=T(k+1,4);

end

 

3.3 结果分析和讨论：

输出结果：

 

44.2800              0              0              0              0

44.2075   44.2026              0              0              0

44.2070   44.2070   44.2070              0              0

44.2070   44.2070   44.2070   44.2070              0

44.2070   44.2070   44.2070   44.2070   44.2070

最终计算得到椭圆周长为44.2070.

如需要代码请在此链接下下载：https://download.csdn.net/download/weixin_41788456/10854229

如有志同者请发邮箱chenshuai0614@hrbeu.edu.cn联系我！