Witam, potrzebuję pilnie użyć funkcji obliczania całki metodą Romberga w Octave lub Matlab. Znalazłem w sieci kilka takowych funkcji, jednakże nie mam o tym zielonego pojęcia, a fukcje pewnie trzeba jakoś wywołać w głównym programie. Byłby bardzo wdzięczny za (nie oszukujmy się) gotowy kod.
Oto funkcje znalezione przeze mnie:
function romberg(f,a,b,n)
% Programista Janusz Chojnacki
%Funkcja romberg oblicza całki oznaczone metodą Romberga
fprintf('\n')
disp(' Tablica Romberga')
disp('_________________________________________________________')
disp(' i h Ri,1 Ri,2 Ri,3 ... ')
disp('_________________________________________________________')
h=b-a;
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
fprintf('%2.0f %8.4f %12.4f\n',1,h,R(1,1));
m=1;
for i=1:n-1
h=h/2;
S=0;
for j=1:m
x=a+h*(2*j-1);
S=S+feval(f,x);
end
4
R(i+1,1)=R(i,1)/2+h*S;
fprintf('%2.0f %8.4f %12.4f',i+1,h,R(i+1,1));
m=2*m;
for k=1:i
R(i+1,k+1)=R(i+1,k)+(R(i+1,k)-R(i,k))/(4^k-1);
fprintf('%12.4f',R(i+1,k+1));
end
fprintf('\n');
end
function Rnum = RombergInt(f,a,b,m)
% Integration algorithm based on Romberg extrapolation
% f - string input for function y = f(x) (e.g. f = 'x.^6')
% a - lower limit of integration
% b - upper limit of integration
% m - maximal number of Romberg iterations
% Rnum - row-vector of numerical approximations for the integral of f(x) between x = a and x = b:
% the entry with index k in Rnum corresponds to the approximation of order O(h^(2k))
R = ones(m,m); % the matrix for Romberg approximations
hmin = (b-a)/2^(m-1); % the minimal step size
for k = 1 : m
h = 2^(k-1)*hmin; % the step size for the k-th row of the Romberg matrix
x = a : h : b;
f1 = eval(f);
k1 = length(f1);
R(k,1) = 0.5*h*(f1(1) + 2*sum(f1(2:k1-1)) + f1(k1)); % trapezoidal rule for the first column of the Romberg matrix
end
for k = 2 : m % compute the k-th column of the Romberg matrix
for kk = 1 : (m-k+1) % the matrix of Romberg approximations is triangular!
R(kk,k) = R(kk,k-1)+(R(kk,k-1) - R(kk+1,k-1))/(4^(k-1)-1); % see the Romberg integration algorithm
end
end
% define the vector Rnum for numerical approximations
Rnum = R(1,:);
function R = romberg(f,a,b,q)
% ROMBERG Integration using Romberg's method
% Usage: R = ROMBERG(FNAME,A,B,Q)
% R is the Romberg table with Q+2 rows (default Q=2)
% for the integral of F=FNAME(X) from A to B
if nargin<4, q = 2; end
R = NaN(q+2,q+2);
h = b-a; m = 1;
R(1,1) = h*(feval(f,a)+feval(f,b))/2;
for i = 2:q+2
h = h/2; m = 2*m;
R(i,1) = R(i-1,1)/2 + h*sum(feval(f,a+h*(1:2:m)));
for j = 1:i-1
R(i,j+1) = R(i,j) + (R(i,j) - R(i-1,j))/(2^(2*j)-1);
end
end