Could you check (if statment block) in my MATLAB code please?

[Avan]

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);
figure(1)
plot(t,p)
xlabel('Time')
ylabel('p')

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 )) ) ) ;

I1 = integral( f, 0,1,'ArrayValued',true)./1
I2 = integral( f, 0,2,'ArrayValued',true)./2
I3 = integral( f, 0,3,'ArrayValued',true)./3
I4 = integral( f, 0,4,'ArrayValued',true)./4

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

if (v < 1e-17)
f = 0
else
f = I4
end

I=[I1,I2,I3,I4] ;
T=[1,2,3,4] ;
figure(2)
plot(T,I)
title('The Fisher Information with time')
xlabel('Time')
ylabel('Fisher Information')

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))
xlabel('Time')
ylabel('The integrand')
1;
% function dpdt = logisticOscilnumerical(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end

The problem as I understood from the code is with the term ( sin(omega*t)) in the denominator ,,I tried to put if statement to avoid singularity but it does not work,, any idea to solve this problem please??

 

[jtbell]

What language is this?

I tried to put if statement to avoid singularity but it does not work

If you show us exactly what you tried, and what you mean by "it does not work" someone can probably suggest a (perhaps small) change to make it work properly.

By the way, your code will be easier to read if you put it between [ code ] and [ \code ] tags. Remove the spaces that are between the brackets [ ]. I had to put them there to make the tags visible here.

 

[Avan]

What language is this?
If you show us exactly what you tried, and what you mean by "it does not work" someone can probably suggest a (perhaps small) change to make it work properly.

By the way, your code will be easier to read if you put it between [ code ] and [ \code ] tags. Remove the spaces that are between the brackets [ ]. I had to put them there to make the tags visible here.

Many thanks.
I am using MATLAB to write my code.
I added the if statement to the code.
sorry about your last suggest, I could not understand what do you mean exactely.

Regards

 

[jtbell]

I am using MATLAB to write my code.

It helps people reading the forum if you put the language name in the thread title. I've edited it, as you can probably see now.

I added the if statement to the code.

Yes, but where? Show us! :)

sorry about your last suggest, I could not understand what do you mean exactely.

I've made a picture to illustrate how to use the tags. I used a C++ program, but the language doesn't matter. When you do this, it will preserve any spacing or indentation that you use to format your code to make it easier to read, instead of collapsing everything to the left.

gives you this:

#include <iostream>

using std::cout;
using std::endl;

int main ()
{
    cout << "Hello, world!" << endl;
    return 0;
}

instead of this:

#include <iostream>

using std::cout;
using std::endl;

int main ()
{
cout << "Hello, world!" << endl;
return 0;
}

 

[Avan]

It helps people reading the forum if you put the language name in the thread title. I've edited it, as you can probably see now.
Yes, but where? Show us! :)

Sorry,I will make it more clear
I've made a picture to illustrate how to use the tags. I used a C++ program, but the language doesn't matter. When you do this, it will preserve any spacing or indentation that you use to format your code to make it easier to read, instead of collapsing everything to the left.

View attachment 76037

gives you this:

#include <iostream>

using std::cout;
using std::endl;

int main ()
{
    cout << "Hello, world!" << endl;
    return 0;
}

instead of this:

#include <iostream>

using std::cout;
using std::endl;

int main ()
{
cout << "Hello, world!" << endl;
return 0;
}

Thanks a lot for your explanation and sorry for that but this is the first time I am posting here.

 

[Nugatory]

Thanks a lot for your explanation and sorry for that but this is the first time I am posting here.

That's cool - now you know how to do it.

You will get more and more helpful answers if you can tell us more about the problem that you're trying to solve. What is this code supposed to do? What is the physical significance of that quantity that you're taking the sine of? If you didn't have a computer to help you with the calculation, were doing it the hard way with pencil and paper, what forumlas would you use if that quantity were zero?

 

[FactChecker]

Which line aborts? Is it in the definition of f? Can you do your analysis from that point on without starting at time 0?
You might try this:
[t_temp,p_temp] = ode45(@logisticOscilnumerical,tspan,p0,[],omega,k,N0);
t=t_temp(2:end);
p=p_temp(2:end);

P.S. There is a separate section for math and science languages where MATLAB experts may help. I don't know how you can move this question there.

 

[Avan]

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);
figure(1)
plot(t,p)
xlabel('Time')
ylabel('p')

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 )) ) ) ;

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

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

if (v < 1e-17)
f = 0
else
f = I4
end

I=[I1,I2,I3,I4] ;
T=[1,2,3,4] ;
figure(2)
plot(T,I)
title('The Fisher Information with time')
xlabel('Time')
ylabel('Fisher Information')

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))
xlabel('Time')
ylabel('The integrand')
1;
% function dpdt = logisticOscilnumerical(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end

The problem as I understood from the code is with the term ( sin(omega*t)) in the denominator ,,I tried to put if statement to avoid singularity but it does not work,, any idea to solve this problem please??

** Note : I do not have any problem with ODE45 , The problem is with f and exactly with the term (sin(omega*t)) in the denominator ,, When I am trying to find I1,I2,I3 I do not have any problem but I could not compute I4 because there is a singularity at (pi = 3.14 ) .
I wrote a ( IF statement ) to avoid integration at this specific point but I think the block of (IF statement ) is not correct because it did not make any change to my result.

 

[Avan]

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);
figure(1)
plot(t,p)
xlabel('Time')
ylabel('p')

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 )) ) ) ;

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

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

if (v < 1e-17)
f = 0
else
f = I4
end

I=[I1,I2,I3,I4] ;
T=[1,2,3,4] ;
figure(2)
plot(T,I)
title('The Fisher Information with time')
xlabel('Time')
ylabel('Fisher Information')

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))
xlabel('Time')
ylabel('The integrand')

1;
% function dpdt = logisticOscilnumerical(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end

The problem as I understood from the code is with the term ( sin(omega*t)) in the denominator ,,I tried to put if statement to avoid singularity but it does not work,, any idea to solve this problem please??

** Note : I do not have any problem with ODE45 , The problem is with f and exactly with the term (sin(omega*t)) in the denominator ,, When I am trying to find I1,I2,I3 I do not have any problem but I could not compute I4 because there is a singularity at (pi = 3.14 ) .
I wrote a ( IF statement ) to avoid integration at this specific point but I think the block of (IF statement ) is not correct because it did not make any change to my result.

 

[DrClaude]

Removing an infinity by defining it away is not a good way to deal with it.

Why do you think that having an if block in the middle of the program will mean that the program is going to refer to it when calling the function f? The if would need to be inside the function.

 

[Avan]

Removing an infinity by defining it away is not a good way to deal with it.

So what do you think ,What is the correct way to deal with this kind of problem please?

Why do you think that having an if block in the middle of the program will mean that the program is going to refer to it when calling the function f? The if would need to be inside the function.

Many thanks.
I tried to put it inside but it still not working ,, I still got a singularity at x = 3.14159.

This is what I did :

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-17)
f = 0
else
I = integral( f, 0.01,4,'ArrayValued',true)./4

end
figure(2)
plot(t,I)1;

 

[DrClaude]

That's not inside the function. The function is defined in a single line in your code. You will need to expand the function in its own MATLAB file.

 

[Avan]

That's not inside the function. The function is defined in a single line in your code. You will need to expand the function in its own MATLAB file.

Could you show me how to do that please?? I do not have a lot of experience with MATLAB.

 

[FactChecker]

** Note : I do not have any problem with ODE45 , The problem is with f and exactly with the term (sin(omega*t)) in the denominator ,, When I am trying to find I1,I2,I3 I do not have any problem but I could not compute I4 because there is a singularity at (pi = 3.14 ) .
I wrote a ( IF statement ) to avoid integration at this specific point but I think the block of (IF statement ) is not correct because it did not make any change to my result.

Sorry. I didn't see the other problems with I4.
Try something like
f=I4;
f(abs(v)<1e-17) = 0;
This will set any value f(I)=0 where abs(v(I)) < 1e-17. You might want to use something like 0.0001 instead of 1e-17 just to be safe.

 

[DrClaude]

Could you show me how to do that please?? I do not have a lot of experience with MATLAB.

http://www.mathworks.com/help/matlab/ref/function.html

 

[Avan]

http://www.mathworks.com/help/matlab/ref/function.html

You mean like this please:

function y = probdenfun(t,p,A,N0,omega,k)
N0=1;
omega=1;
p0=0.1;
k=10;
t=0:0.1:4;
y = ( ( 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-17)
y = 0
else
f = y
end
end

 

[DrClaude]

I don't understand why you set t inside the function. Shouldn't it be a parameter of the function?

You should calculate v first, then check its value, and then calculate y. Otherwise, you are still doing a division by zero, which can cause you some problems. Also, there should be no variable f in the function, as you are not using it.

 

[Avan]

Sorry. I didn't see the other problems with I4.
Try something like
f=I4;
f(abs(v)<1e-17) = 0;
This will set any value f(I)=0 where abs(v(I)) < 1e-17. You might want to use something like 0.0001 instead of 1e-17 just to be safe.

Many thanks but unfortunately, it does not work.

 

[Avan]

I don't understand why you set t inside the function. Shouldn't it be a parameter of the function?

You should calculate v first, then check its value, and then calculate y. Otherwise, you are still doing a division by zero, which can cause you some problems. Also, there should be no variable f in the function, as you are not using it.

Sorry,I could not understand what do you mean. Can you show me where is the error that I am doing please
Thanks

 

[DrClaude]

function y = probdenfun(t,p,A,N0,omega,k)
   N0=1;  % ins't N0 a parameter of the function?
   omega=1; % ins't omega a parameter of the function?
   p0=0.1; % p0 is not used
   k=10; % ins't k a parameter of the function?
   t=0:0.1:4; % ins't t a parameter of the function?

% don't calculate y yet, do v first
   y = ( ( 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-17)
      y = 0
   else
% calculate y here
      f = y  % f is not used 
   end
end

 

[Avan]

function y = probdenfun(t,p,A,N0,omega,k)
   N0=1;  % ins't N0 a parameter of the function?
   omega=1; % ins't omega a parameter of the function?
   p0=0.1; % p0 is not used
   k=10; % ins't k a parameter of the function?
   t=0:0.1:4; % ins't t a parameter of the function?

% don't calculate y yet, do v first
   y = ( ( 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-17)
      y = 0
   else
% calculate y here
      f = y  % f is not used
   end
end

Now it is correct please?

function y = probdenfun(t,p,A,N0,omega,k)

v = ( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k).^2 ))

if (v < 1e-17)
y = 0
else
y = ( ( A.*( ( N0.* (sin(omega.*t)).^2 .*(1-(2.*p./k))+(omega.*cos(omega.*t) ) ).^2 ) ./( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k)).^2 )) );
end
end

 

[jtbell]

P.S. There is a separate section for math and science languages where MATLAB experts may help. I don't know how you can move this question there.

Oops, I should have done that in the first place. I'll do it now...

[added] I found a second thread on this topic in the "Math and Science Textbooks" forum, which is not the same as the "Math and Science Software" forum. The two threads have now been merged into a single thread in the software forum.

 

[FactChecker]

Many thanks but unfortunately, it does not work.

When I ran a small test, it did exactly what your 'if' test was trying to do. If you are not willing to indicate more than "it doesn't work" than I can't help. Is there an error message? What line does the message indicate? Please be more informative. That is the least you can do.

 

[Avan]

When I ran a small test, it did exactly what your 'if' test was trying to do. If you are not willing to indicate more than "it doesn't work" than I can't help. Is there an error message? What line does the message indicate? Please be more informative. That is the least you can do.

I am sorry for that and many thanks for your help.

What I mean by it does not work that I still received a warning message about singularity at 3.14 after including if statement. But if it is work with you so could you show me how you used it please?

 

[Avan]

Now it is correct please?

function y = probdenfun(t,p,A,N0,omega,k)

v = ( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k).^2 ))

if (v < 1e-17)
y = 0
else
y = ( ( A.*( ( N0.* (sin(omega.*t)).^2 .*(1-(2.*p./k))+(omega.*cos(omega.*t) ) ).^2 ) ./( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k)).^2 )) );
end
end

Now, I have two MATLAB file as follow:

The first one is :

function dpdt = logisticOscilnumerical(t,p,omega,k,N0)
dpdt = N0*sin(omega*t)*p*(1-p/k);
end

The second one is :

function y = probdenfun(t,p,A,N0,omega,k)
v = ( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k).^2 ))
if (v < 0.0001)
y = 0
else
y = ( ( A.*( ( N0.* (sin(omega.*t)).^2 .*(1-(2.*p./k))+(omega.*cos(omega.*t) ) ).^2 ) ./( (N0.^2).*((sin(omega.*t)).^4).*((p-p.^2./k)).^2 )) );
end
end

Now how could I call both of them in my code please?

First I need to solve the system numerically using ODE45.
Second , Finding the Integral for the y function which is a PDF for the system.
Because I got a singularity so I tried to avoid it using if statement.

Regards

 

[FactChecker]

An 'if' condition doesn't work with an entire array. Here is an example of how to find 0 values in one array and set the corresponding values of another array:

A = [ 1 2 3 0 5] % array with one problem value at A(4)
B = [ 6 7 8 9 10] % another series of numbers matched to A
B(abs(A)<0.1) = 0 % for any entry of A with abs(A(i)) < 0.1, set B(i)=0

The results are:

A =
  1  2  3  0  5
B =
  6  7  8  9  10
B =
  6  7  8  0  10

I am suspicious that the problem occurs before your 'if' code. I don't see how an integration can be done through the singularity.
I recommend that you check what line of MATLAB code your error message is coming from and look at that in detail.

 

[Avan]

An 'if' condition doesn't work with an entire array. Here is an example of how to find 0 values in one array and set the corresponding values of another array:

A = [ 1 2 3 0 5] % array with one problem value at A(4)
B = [ 6 7 8 9 10] % another series of numbers matched to A
B(abs(A)<0.1) = 0 % for any entry of A with abs(A(i)) < 0.1, set B(i)=0

The results are:

A =
  1  2  3  0  5
B =
  6  7  8  9  10
B =
  6  7  8  0  10

I am suspicious that the problem occurs before your 'if' code. I don't see how an integration can be done through the singularity.
I recommend that you check what line of MATLAB code your error message is coming from and look at that in detail.

Do you mean using like this block with my code?

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));
a = ones(401,1);
if a(abs(v)<0.0001) == 0

The problem is finding I4 .

 

[FactChecker]

The problem is finding I4 .

That is what has me confused. In your original post, you were trying to correct f after I4 had already been calculated. I don't see how that would help the I4 calculation. That is why I asked what MATLAB line is aborting. (I don't have MATLAB at home so I can't run your code.)

As far as I know, you can not put a comparison of vector elements in the 'if' condition. "if v < 0.000001 " does not work element-by-element.
Here is one way to avoid that. Replace line 15 where f is defined with this code:

my_lower_limit = 0.0000001;
my_sin = sin(omega.*t);

% keep positive sin away from 0 and replace with a small positive number
my_sin( (0 <= my_sin) && (my_sin < my_lower_limit) ) = my_lower_limit;
% keep negative sin away from 0 and replace with a small negative number
my_sin( (my_sin < 0) && (my_sin > -my_lower_limit) ) = -my_lower_limit;

f = @(t) ( ( A.*( ( N0.* (my_sin).^2 .*(1-(2.*P(t)./k))+(omega.*cos(omega.*t) ) ).^2 ) ./( (N0.^2).*((my_sin).^4).*((P(t)-(P(t).^2./k)).^2 )) ) ) ;

Set the value of my_lower_limit small enough for the accuracy you need.

Another way that may be simpler is to loop through the individual elements of v and apply your logic to each one.
Replace the definition of f in line 15 with this:

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

% get number of rows in the column vector v
[num_rows,num_columns] = size(v);

% loop through the values of v and define f_vector appropriately
for i = 1:num_rows
   if abs(v(i)) < 1e-17
      f_vector(i) = 0;
   else
      f_vector(i) = ( ( A*( ( N0* (sin(omega*t(i)))^2 *(1-(2*P(t(i))/k))+(omega*cos(omega*t(i)) ) )^2 )/v(i) ) );
   end
end

% Define f.  ? I'm not familiar enough with this to know if it is correct
f = @(t) f_vector;

Since I don't have MATLAB here, I can't tell if I got the syntax right. You may have to correct parts.

 

[Avan]

That is what has me confused. In your original post, you were trying to correct f after I4 had already been calculated. I don't see how that would help the I4 calculation. That is why I asked what MATLAB line is aborting. (I don't have MATLAB at home so I can't run your code.)

As far as I know, you can not put a comparison of vector elements in the 'if' condition. "if v < 0.000001 " does not work element-by-element.
Here is one way to avoid that. Replace line 15 where f is defined with this code:

my_lower_limit = 0.0000001;
my_sin = sin(omega.*t);

% keep positive sin away from 0 and replace with a small positive number
my_sin( (0 <= my_sin) && (my_sin < my_lower_limit) ) = my_lower_limit;
% keep negative sin away from 0 and replace with a small negative number
my_sin( (my_sin < 0) && (my_sin > -my_lower_limit) ) = -my_lower_limit;

f = @(t) ( ( A.*( ( N0.* (my_sin).^2 .*(1-(2.*P(t)./k))+(omega.*cos(omega.*t) ) ).^2 ) ./( (N0.^2).*((my_sin).^4).*((P(t)-(P(t).^2./k)).^2 )) ) ) ;

Set the value of my_lower_limit small enough for the accuracy you need.

Another way that may be simpler is to loop through the individual elements of v and apply your logic to each one.
Replace the definition of f in line 15 with this:

v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

% get number of rows in the column vector v
[num_rows,num_columns] = size(v);

% loop through the values of v and define f_vector appropriately
for i = 1:num_rows
   if abs(v(i)) < 1e-17
      f_vector(i) = 0;
   else
      f_vector(i) = ( ( A*( ( N0* (sin(omega*t(i)))^2 *(1-(2*P(t(i))/k))+(omega*cos(omega*t(i)) ) )^2 )/v(i) ) );
   end
end

% Define f.  ? I'm not familiar enough with this to know if it is correct
f = @(t) f_vector;

Since I don't have MATLAB here, I can't tell if I got the syntax right. You may have to correct parts.

Many thanks for your help.
The second suggestion is working but I am getting something different. I will explain what I mean as follow, So what do you think why there is a different please?
1- For my original code when I tried to plot (T,I) where T is from 0.01 to 3 and without a singularity problem, the code and plot is as following :

function RunlogOscil
omega=1;
N0=1;
k = 10;
A = 1;
p0 = 0.1;
tspan=(0:0.01:3);

[t,p] = ode45(@logOscilnumeric,tspan,p0,[],omega,k,N0);

figure(1)
plot(t,p)

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 )) ) ;
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
I=[I1,I2,I3];
T=[1,2,3];

figure(2)
plot(T,I)

P = @(T) interp1(t,p,T);
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 )) )
f
igure(3)
plot(tspan,g(tspan))

1;
% function dpdt = logOscil(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p*(1-p/k);
% end

 

[Avan]

2- After I added your suggestion to my code to avoid a singularity when I am integration from 0.01 to 4, I got different plot as follow :

function RunlogOscilnumerfishfixedn0omega

omega=1;
N0=1;
k = 10;
A = 1;
p0 = 0.1;
tspan=(0:0.01:4);

[t,p] = ode45(@logOscilnumeric,tspan,p0,[],omega,k,N0);

figure(1)
plot(t,p)

P = @(T) interp1(t,p,T);
v = ( (N0.^2).*((sin(omega.*t)).^4).*((P(t)-(P(t).^2./k)).^2 ));

% get number of rows in the column vector v
[num_rows,num_columns] = size(v);
% loop through the values of v and define f_vector appropriately
for i = 1:num_rows
if abs(v(i)) < 1e-17
f_vector(i) = 0;
else
f_vector(i) = ( ( A*( ( N0* (sin(omega*t(i)))^2 *(1-(2*P(t(i))/k))+(omega*cos(omega*t(i)) ) )^2 )/v(i) ) );
end
end

f = @(t) f_vector(i);

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)

P = @(T) interp1(t,p,T);
g = @(t) ( ( A*( ( N0* (sin(omega*t(i)))^2 *(1-(2*P(t(i))/k))+(omega*cos(omega*t(i)) ) )^2 )/v(i) ) );

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