CODEWITHSUNDEEP
matlabwolfram-mathematica
WnGatRC456
WnGatRC456

Reputation: 337

Evaluate indefinite integral numerically in matlab/mathematica that it cannot do symbolically

I am trying to calculate the integral of a function in Matlab and Mathematica that the software cannot do symbolically.

Here is my MatLab code so far, but I understand it may not be very helpful as is.

f = @(t) asin(0.5*sin(t));
a = @(t) sin(t);
F = int(f,t)   % Matlab can't do this
F = 
int(asin(sin(t)/2), t)
A = int(a,t)   % This works
A =
-cos(t)

dt = 1/(N-1); % some small number
for i=1:N
    F(i) = integral(f,(i-1)*dt,i*dt);
    A(i) = integral(a,(i-1)*dt,i*dt);
end

Both of the calculations in the for loop give a rough approximation of f or a not their integrals after multiplying by dt.

On the math stack-exchange I found a question that derives a finite difference like method for the integral at a point. However, when I did the calculation in Matlab it output a scaled down version of f which was evident after plotting (see above for what I mean by scaled down). I think that's because for smaller intervals the integral basically approximates the function to varying degrees of accuracy (again see above).

I am trying to get either a symbolic equation for the integral, or an approximation of the integral of the function at each location.

So my question is then if I have a function f that MatLab and Mathematica cannot easily take the integral of

  1. can I approximate the integral directly with an integral calculator besides the default ones? (int,integral,trapz)

or

  1. can I approximate the function with finite differences first and then evaluate the integral symbolically?

Upvotes: 1

Views: 692

Answers (2)

WnGatRC456
WnGatRC456

Reputation: 337

The accepted answer in general is by far the best method I would say but if certain restrictions on your functions are allowable then there is a second method.

For two functions f and g see below

T = 1;  % Period
NT = 1;  % Number of periods
dt = 0.01; % time interval
time = 0:dt:NT*T;  % time

syms t
x = K*sin(2*pi*t+B);   % edit as appropriate

% f = A/tanh(K)*tanh(K*sin(2*pi*t+p))
% g = A/asin(K)*asin(K*sin(2*pi*t+p))

formulas found here

f = A1/tanh(K1)*(2^(2*1)-1)*2^(2*1)*bernoulli(2*1)/factorial(2*1)*x^(2*1-1);
% |K1|<pi/2
g = A2/asin(K2)*factorial(2*0)/(2^(2*0)*factorial(0)^2*(2*0+1))*x^(2*0+1);
% |K2|<1

there are no such limitations in the accepted answer

N = 60;
for k=2:N
    a1 = (2^(2*k)-1)*2^(2*k)*bernoulli(2*k)/factorial(2*k);
    f = f + A1/tanh(K1)*a1*x^(2*k-1);

    a2 = factorial(2*k)/(2^(2*k)*factorial(k)^2*(2*k+1));
    g = g + A2/asin(K2)*a*x^(2*k+1);
end

MATLAB can calculate sin^n(t) for n being an integer.

F = int(f,t);
phi = double(subs(F,t,time));

G = int(g,t);
psi = double(subs(G,t,time));

Upvotes: 0

PilouPili
PilouPili

Reputation: 2699

Your code is nearly fine it's just that

for i=1:N
    F(i) = integral(f,0,i*dt);
end

You could also do

F(1)=integral(f,0,dt)
for i=2:N
    F(i) = F(i-1)+integral(f,(i-1)*dt,i*dt);
end

Second option is surely more efficient

Because the primitive is really F(x)=int(f(x), 0, x) (0 defines a certain constant ) and for sufficiently small dx you have shown that f(x)=int(f(x), x,x+dx)/dx i. You have proven that MATLAB intégral function does its job.

For example let's take enter image description here= enter image description here the function above will compute enter image description here if you wish to compute enter image description here just replace 0 above by the constant a you like.

now enter image description here and so you should get F containing a discretization of enter image description here

Upvotes: 1

Related Questions