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

1 次查看（过去 30 天）
Avan Al-Saffar
Avan Al-Saffar 2014-10-24
评论： Stephen23 2015-8-28
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

请先登录，再回答此问题。

回答（4 个）

Torsten
Torsten 2014-10-24
Use MATLAB's dsolve to solve your ODE analytically instead of ODE45.
I guess this will resolve your integration problems.
Best wishes
Torsten.
Torsten
Torsten 2014-10-27
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 2014-10-27
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 个评论
显示 6更早的评论
Torsten
Torsten 2014-12-1
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 2014-12-1
编辑：Avan Al-Saffar 2014-12-1
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

请先登录，再进行评论。

Torsten
Torsten 2014-12-1
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 个评论
Avan Al-Saffar
Avan Al-Saffar 2014-12-1
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by