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

Spruce budworm diff.eq,parameter variation using matlab functions

  1. Dec 13, 2012 #1
    I have come across a matlab code for solving the spruce budworm differential equation.But,I would like to solve the same differential equation for a range of parameters(r=0:5,q=0:10).I am having problems trying to define the array of matrices for each loop.Please look at the code below:

    % This program solves the differential equation
    % du/dt = f(t,u) = ru(1-u/q) + u^2/(1 + u^2)

    function spruceco

    tspan = [0; 100]; %This command defines the time interval of interest.
    u0 = [.1];
    for q=0:10 %I would like to add these here and then pass the values of q,r to the function f%
    for r=0:5

    [t,u] = ode15s(@f,tspan,u0)
    % This command tells MATLAB to use the differential equation solver called ode15s to numerically compute the solution of the differential equation defined by the function f, for time interval tspan, and initial condition u0. The right-hand side of the equation tells MATLAB to store this output in vectors t (for the time points) and u (for the population density).

    plot(t,u(:,1)); %This tells MATLAB to plot the solution.
    hold on
    % --------------------------------------------------------------------------

    function dudt = f(t,u) %This commands defines the function du/dt = f(t,u).
    q=10;r=2;%I would like to solve the differential equation for q=0:10,r=0:5;How?Procedure?
    dudt = [r*u(1)*(1-u(1)/q) - u(1)^2/(1+u(1)^2)]; %This command inputs the left-hand side of the spruce budworm
    %differential equtaion.
  2. jcsd
  3. Dec 17, 2012 #2
    I'm not a Matlab user, but I think the technique of adding r and q to the differential equation might work. Essentially, initialize u0 to be a vector with values 0.1 (your current value), r and q. You then define the differential equation f(t,u) to return it's current value as the first element of a 3-vector and 0 for the other 2 elements (r and q being constant throughout each integration). You can then loop through r and q values, something like the following (or however Matlab constructs vectors).

    f(t,u) = [u(2).u(1)(1-u(1)/u(3)) + u(3)^2/(1 + u(3)^2) 0 0]

    for r = ....
    for q = ....
    u0 = [0.1 r q]
    [t,u] = ode15s(@f,tspan,u0)

    The image shows something similar in Mathcad
    http://https://www.physicsforums.com/attachment.php?attachmentid=54011&stc=1&d=1355766195 [Broken]
    You also might want to check that the equation is valid for your desired r and q ranges ... (eg, what the result if q = 0 ?)

    Attached Files:

    Last edited by a moderator: May 6, 2017
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook