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

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

[Stephen23]

@Avan Al-Saffar: please do not put empty lines between every line of code.

Answer 1

[Torsten]

Use MATLAB's dsolve to solve your ODE analytically instead of ODE45.

I guess this will resolve your integration problems.

Best wishes

Torsten.

[Avan Al-Saffar]

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.

Answer 2

[Torsten]

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.

[Avan Al-Saffar]

Dear Torsten Thanks again.Could you show me how please? Because, I am not sure that I understood what you mean?

Regards

Avan

Answer 3

[Torsten]

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.

[Torsten]

What I mean is that

P(t)=1/(1/P0+1/k*(1-exp(-N0/omega*(1-cos(omega*t)))))

is the analytical solution for your differential equation

dP/dt=N0*sin(omega*t)*P*(1-P/k) , P(0)=P0.

So you can work with this explicit solution instead of using the numerical obe obtained from ODE45.

Best wishes

Torsten.

[Avan Al-Saffar]

Many thanks, I tried to work with your suggestion but I still have the same problem.Any other advice please? Actually, I want to ask you if could I depend on Analytical results instead of numerical results ?? I am asking you because this system is a simple system which could I solved analytically what is about the complex systems??

Regards

[Torsten]

Concerning other advice:

I don't know where your function f comes from. So I can't give you advice on this.

Concerning analytical versus numerical solutions:

For better precision, you should use analytical solutions as long as you can.

Only switch to numerical solutions if analytical solutions are no longer available or too complicated to be evaluated efficiently.

Best wishes

Torsten.

[Avan Al-Saffar]

I will give you more details :

u = dx/dt = N0*sin(omega*t)*(p-p^2/k)

What I want to find is this form :

I = (1/T)*integral(f,0,T)

f = (1/u^4)*(du/dt)^2

du/dt = N0 * (p-p^2/k)*(omega*cos(omega*t) + N0 * sin(omega*t)^2*(1-2*p/k)

Regards

[Torsten]

I thought that maybe your f looks like

1/p^4*(dp/dt)^2

But if it's formed with the derivatives

(1/(dp/dt)^4)*(d^2p/dt^2)^2

- that's bad luck!

Best wishes

Torsten.

[Avan Al-Saffar]

So this means that I Should ignore all your comments?

Regards

[Torsten]

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]

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

Answer 4

[Torsten]

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.

[Avan Al-Saffar]

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