Could anyone help me with tolerances erroe that I got when I am trying to implement an integration please??

1 view (last 30 days)
Avan Al-Saffar
Avan Al-Saffar on 24 Oct 2014
Commented: Stephen23 on 28 Aug 2015
function RunLogisticOscilFisher
omega=1;
k=10;
N0=1;
A=1;
p0=.1;
tspan=(0:0.01:100);
% Finding the numerical solution for the function using ode45 solver
[t,p]=ode45(@logisticOscilfisher,tspan,p0,[],N0,k,omega);
% Plotting the function with time
figure(1)
plot(t,p)
P = @(T) interp1(t,p,T);
% Finding the integral to get the Fisher Information
f = @(T) ( A*(((N0*sin(omega*T).^2.*(1-2*P(T)./k))+(omega.*cos(omega*T) ) ).^2)./(N0.^2*sin(omega*T).^4.*(P(T)-P(T).^2./k).^2) )
I1=integral( f, 1,10)
I2=integral( f, 1,40,'ArrayValued',true)
I3=integral( f, 1,60,'ArrayValued',true)
I4=integral(f,1,80,'ArrayValued',true)
I5=integral(f,1,100,'ArrayValued',true)
I=[I1./20 I2./40 I3./60 I4./80 I5./100]
R=[20 40 60 80 100];
%Plotting the Fisher Information
figure(2)
plot(R,I);
1;
% function dpdt = logisticOscilfisher(t,p,N0,k,omega)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end

Answers (4)

Torsten
Torsten on 24 Oct 2014
Use MATLAB's dsolve to solve your ODE analytically instead of ODE45.
I guess this will resolve your integration problems.
Best wishes
Torsten.
  1 Comment
Avan Al-Saffar
Avan Al-Saffar on 26 Oct 2014
Many thanks for your answer but actually, I do not have a problem with ODE45. I am getting an error with the integration which is related to finding I.

Sign in to comment.


Torsten
Torsten on 27 Oct 2014
Of course I'm not sure, but I guess the problem with the tolerances stems from the interpolation of the solution obtained from ODE45 within your function f to be integrated. Thus having an explicit expression for P will help for the integration. MATLAB's dsolve will give you this explicit expression.
Best wishes
Torsten.
  1 Comment
Avan Al-Saffar
Avan Al-Saffar on 27 Oct 2014
Dear Torsten Thanks again.Could you show me how please? Because, I am not sure that I understood what you mean?
Regards
Avan

Sign in to comment.


Torsten
Torsten on 27 Oct 2014
Try
P=@(T)(1./(1/p0+1/k*(1-exp(-N0/omega*(1-cos(omega*T))))));
instead of
P = @(T) interp1(t,p,T);
in your code above.
Best wishes
Torsten.
  8 Comments
Avan Al-Saffar
Avan Al-Saffar on 1 Dec 2014
thank you to be patient. But I am getting a result for I unless at 0,pi,2*pi,...
I am really confused now.
Please I will ask you just to be more clear, I have a system and I am using ode45 to solve it so do you think there is a problem here?
Secondly, I am finding an integral for f which is as mentioned in my code,, what is about this step please?
Regards

Sign in to comment.


Torsten
Torsten on 1 Dec 2014
There is no problem solving the ODE
dp/dt = N0*sin(omega*t)*p*(1-p/k)
using ODE45, I guess.
You can easily check this by deleting everything below the line
plot(t,p)
in your code above.
The problem is the function f you try to integrate from 1 to 10, from 1 to 40 etc.
It has singularities at points pi,2*pi,3*pi,... (the denominator is zero) and I1, I2, I3,... do not exist (are infinity in this case).
I don't know what you mean by "I am getting a result for I unless at 0,pi,2*pi,...". Could you clarify ?
Best wishes
Torsten.
  1 Comment
Avan Al-Saffar
Avan Al-Saffar on 1 Dec 2014
Dear Torsten
Many thanks.
You are definitely right( in another thread that I submitted I tried to find the integration from 1 to 2 and from 1 to 3 so I am getting values, this is what I mean.
As I know, for every problem there is a solution and this is what I tried to do with putting ( if statement ) with this code in another thread but unfortunately, I do not know how to construct it correctly so could you offer some help with the another thread that I submitted please?
Regards

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!