Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Dec 1, 2014 #1
    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. jcsd
  3. Dec 1, 2014 #2

    jtbell

    User Avatar

    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
  4. Dec 1, 2014 #3
    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
     
  5. Dec 1, 2014 #4

    jtbell

    User Avatar

    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.

    hello.gif

    gives you this:

    Code (Text):

    #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;
    }
     
  6. Dec 1, 2014 #5
    Thanks a lot for your explanation and sorry for that but this is the first time I am posting here.
     
  7. Dec 1, 2014 #6

    Nugatory

    User Avatar

    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?
     
  8. Dec 1, 2014 #7

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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
  9. Dec 2, 2014 #8
    ** 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 #9
    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.
     
  11. Dec 2, 2014 #10

    DrClaude

    User Avatar

    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.
     
  12. Dec 2, 2014 #11
    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;
     
  13. Dec 2, 2014 #12

    DrClaude

    User Avatar

    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.
     
  14. Dec 2, 2014 #13
    Could you show me how to do that please?? I do not have a lot of experience with MATLAB.
     
  15. Dec 2, 2014 #14

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  16. Dec 2, 2014 #15

    DrClaude

    User Avatar

    Staff: Mentor

  17. Dec 2, 2014 #16
    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
     
  18. Dec 2, 2014 #17

    DrClaude

    User Avatar

    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.
     
  19. Dec 2, 2014 #18
    Many thanks but unfortunately, it does not work.
     
  20. Dec 2, 2014 #19
    Sorry,I could not understand what do you mean. Can you show me where is the error that I am doing please
    Thanks
     
  21. Dec 2, 2014 #20

    DrClaude

    User Avatar

    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
     
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Could you check (if statment block) in my MATLAB code please?
Loading...