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

2 Ansichten (letzte 30 Tage)
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

Antworten (4)

Torsten
Torsten am 24 Okt. 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 Kommentar
Avan Al-Saffar
Avan Al-Saffar am 26 Okt. 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.

Melden Sie sich an, um zu kommentieren.


Torsten
Torsten am 27 Okt. 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.

Torsten
Torsten am 27 Okt. 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 Kommentare
Torsten
Torsten am 1 Dez. 2014
It means that if f=1/u^4*(du/dt)^2, I=infinity, independent from whether you supply P analytically or numerically.
Best wishes
Torsten.
Avan Al-Saffar
Avan Al-Saffar am 1 Dez. 2014
Bearbeitet: Avan Al-Saffar am 1 Dez. 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

Melden Sie sich an, um zu kommentieren.


Torsten
Torsten am 1 Dez. 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 Kommentar
Avan Al-Saffar
Avan Al-Saffar am 1 Dez. 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

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by