Why the if statement does not work when I am trying to avoid the singularity in my code ??

1 次查看（过去 30 天）
Avan Al-Saffar
Avan Al-Saffar 2014-11-17
编辑： Avan Al-Saffar 2014-12-3
function RunlogisticOscilnumericalfisherfixedn0omega
omega=1;
N0=1;
k = 10;
A = 1;
p0 = 0.1;
tspan=(0:0.1:4);
[t,p] = ode45(@logisticOscilnumerical,tspan,p0,[],omega,k,N0);
P = @(T) interp1(t,p,T);
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 )) ) ) ;
v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 )); if v <= 1e-100
f = 0
else
f = I4
end
I1 = integral( f, 0.01,1,'ArrayValued',true)./1
I2 = integral( f, 0.01,2,'ArrayValued',true)./2
I3 = integral( f, 0.01,3,'ArrayValued',true)./3
I4 = integral( f, 0.01,4,'ArrayValued',true)./4
I=[I1,I2,I3,I4] ;
T=[1,2,3,4] ;
figure(2)
plot(T,I)
g = @(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 )) ) ) ;
figure(3)
plot(tspan,g(tspan))
1;
% function dpdt = logisticOscilnumerical(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end
Could anyone help me to avoid singularity at x = 3.14159 using if statement please??
  5 个评论
显示 3更早的评论
Avan Al-Saffar
Avan Al-Saffar 2014-11-28
Ok, I know that my problem is with the term(sin(omega*t) but what if I tried to use if statement to avoid integration at those specific values, What do you think? Or do you have another idea?
Regards
Avan
Torsten
Torsten 2014-11-28
The values of a function at a finite number of points do not influence the integral.
Changing the function in a neighborhood of the singularities will give you a finite value for the integrals, but this finite value will depend on the size of the neighborhood you chose.
So I think it is better to reconsider how you came up with the function f and if this f is appropriate for your purpose.
Best wishes
Torsten.

请先登录，再进行评论。

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

回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by