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

1. Dec 1, 2014

### 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??

Last edited: Dec 1, 2014
2. Dec 1, 2014

### Staff: Mentor

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.

Last edited: Dec 1, 2014
3. Dec 1, 2014

### Avan

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

4. Dec 1, 2014

### Staff: Mentor

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! :)

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:

Code (Text):

#include <iostream>

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

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

#include <iostream>

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

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

5. Dec 1, 2014

### Avan

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

6. Dec 1, 2014

### Staff: Mentor

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?

7. Dec 1, 2014

### 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.

Last edited: Dec 1, 2014
8. Dec 2, 2014

### Avan

** 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.

9. Dec 2, 2014

### 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.

10. Dec 2, 2014

### Staff: Mentor

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.

11. Dec 2, 2014

### Avan

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;

12. Dec 2, 2014

### Staff: Mentor

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.

13. Dec 2, 2014

### Avan

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

14. Dec 2, 2014

### FactChecker

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.

15. Dec 2, 2014

### Staff: Mentor

16. Dec 2, 2014

### Avan

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

17. Dec 2, 2014

### Staff: Mentor

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.

18. Dec 2, 2014

### Avan

Many thanks but unfortunately, it does not work.

19. Dec 2, 2014

### Avan

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

20. Dec 2, 2014

### Staff: Mentor

Code (Text):

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

21. Dec 2, 2014

### Avan

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

Last edited: Dec 2, 2014
22. Dec 2, 2014

### Staff: Mentor

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.

Last edited: Dec 2, 2014
23. Dec 2, 2014

### FactChecker

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.

24. Dec 3, 2014

### Avan

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?

25. Dec 3, 2014

### Avan

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