Could anyone help me with tolerances erroe that I got when I am trying to implement an integration please??
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
1 commentaire
Stephen23
le 28 Août 2015
@Avan Al-Saffar: please do not put empty lines between every line of code.
Réponses (4)
Torsten
le 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.
Torsten
le 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.
Torsten
le 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 commentaires
Torsten
le 1 Déc 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.
Torsten
le 1 Déc 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.
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!