romberg_integration.m

来自「Romberg integration algorithm: a method 」· M 代码 · 共 32 行

%%  maryam-mahmoudi    S.N:870874022
%%
clc;clear;close all;
digits(9);
a = input('please enter a   ');
b = input('please enter b   ');
n = input('please enter n   ');
fx= input('please enter fx  ','s');
tol=input('please enter tol   ');
%%
f=inline(fx);
m=1;
h=(b-a);
cls=1;
j=1;
r(1,1)=h*(f(a)+f(b))/2;
while (cls>tol && j<n) || j<4
    j=j+1;
    h=h/2; 
    sum=0; 
    for p=1:m
        sum=sum+f(a+h*(2*p-1));
    end
        r(j,1)=r(j-1,1)/2+h*sum;
        m=2*m;
       for k=2:j
           r(j,k)=r(j,k-1)+(r(j,k-1)-r(j-1,k-1))/(4^k-1);
       end 
cls=abs(r(j-1,j-1)-r(j,j)) ;
end
fprintf  (['the best approximation value of the integral using Romberg integration is    ', int2str(r(j,j))])
vpa(r)