Ode solver varying parameters

Hello,

I have a three differnetial equations with 8 parameters in total of which 1 parameter p7 actually changes with time (p7=a*b^y*c^z - given for representative purpose, but i have a similar equation for that). How do i solve this with ode45? When i tried i got the following error: Unable to perform assignment because the left and right sides have a different number of elements. This is because the p7 has several values wheras other parametrs have only 1. Can someone help me how to solve ode with changing parameters.

CODE:

dx = zeros(3,1);
  p7=a*b^y*c^z;
  dx(1)=((p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))).*x(1));              
  dx(2)=((-((1/p4).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1));
  dx(3)=(((-((1/p6).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1))+(p7.*(p8-x(3))));

Answer 1

[Walter Roberson]

You do not assign multiple values to p7. Instead, each time that p7 needs to change values in a discontinuous way, you need to terminate the ode45, and then restart it from where you left off except with the new p7 value.

If the p7 changes happen at known times, then the way to do this is to call ode45 with tspan reflecting the interval in which p7 is to remain constant. For example,

p7times = [0, 1, 2.5, 5, 5.2, 10];
p7vals = [pi/9, 0.3, exp(-1), sqrt(2), -3];  %one fewer
num_interval = length(p7vals);
x0 = [0 0 0];
t = cell(num_interval,1);
x = cell(num_interval, 1);
for idx = 1 : num_interval
    p7 = p7vals(idx);
    [t{idx},x{idx}] = ode45(@(t,x) try_eqns(t, x, p7), p7times(idx:idx+1), x0);
    x0 = x{idx}(end,:);
end

If discontinuous changes to your p7 depend upon the results of the ode, then you should instead use ode event functions to detect the situation and terminate the integration and then resume with the new p7. There is a useful example of this: see ballode

[Nivedhitha S]

Thank you Mr.Walter, I tried your code.

As you had mentioned, p7 happens at known times and i modified the code suggested by you.

p7times = [0:60:18000]';
p7vals =[%all p7 values goes here: 300 values];
num_interval = length(p7vals);
x0 = [2.07155 20 0.008];
t = cell(num_interval,1);
x = cell(num_interval, 1);
for idx = 1 : num_interval
    p7 = p7vals(idx);
    [t{idx},x{idx}] = ode45(@(t,x) try_eqns(t, x, p7), p7times(idx:idx+1), x0);
    x0 = x{idx}(end,:);
end

But i got one question:

Since t and x are assigned as cell, i get the output as follows: t has 300 cells with 41*1 double type data. ie. every time interval say 0-60 has been divided into 41 intervals and the ode is solved for every t span data. How i do i avoid this? (say, if i want the solver to work only for t=0, t=60, ., so on.) I tried several methods of reassigning the t, changing it to double instead of a cell etc., But couldn't achieve what i wanted. Could you please help me, sir?

[Walter Roberson]

p7times = [0:60:18000]';
p7vals =[%all p7 values goes here: 300 values];
num_interval = length(p7vals);
x0init = [2.07155 20 0.008];
x0 = x0init;
num_vars = length(x0);
t = zeros(num_interval,1);
x = zeros(num_interval, num_vars);
for idx = 1 : num_interval
    p7 = p7vals(idx);
    [this_t,this_x] = ode45(@(t,x) try_eqns(t, x, p7), p7times(idx:idx+1), x0);
    t(idx) = this_t(end);
    x0 = this_x(end,:);
    x(idx,:) = x0;
end

%after the loop, there is information about the _end_ of each interval
%add in the information about the first interval
t = [p7times(1); t];
x = [x0init; x];

[Nivedhitha S]

Thanks Mr.Walter.

To get clear, on what basis are the variable this_t and this_x assigned? Why not have the equations getting solved only fotr the given t spans? I got about 12937 values for these variables. I am not able to follow the actual mechanism which this is done. It would be great if you could help me out.

Thanks a lot.

[Walter Roberson]

At each ode45 call you are passing a vector of two times. When you pass a vector of length 2 then ode45 returns for a series of times according to whatever it feels like internally, and the number of times returned for can vary depending upon integration tolerances (so two calls for slightly different time ranges are not necessarily going to return the same number of outputs.)

When you pass ode45 a vector of 3 or more times then it returns only for those times, except as needed for terminal events.

Either way, ode45 will evaluate wherever it wants, and the question is just which times it reports back about. It is *not* the case that it only evaluates the function at the given times: it has to evaluate at intermediate times to meet tolerances.

If you wanted then you can pass in [t1, (t1+t2)/2, t2] to reduce the number of times it reports back to 3,but you still end up having to extract the (end) output.

You cannot just pass in list of p7 times because you are changing p7 abruptly and ode45 requires that the function be differentiable

Answer 2

[J Chen]

You put the derivative calculation in a file such as fcn.m,

function dx = fcn(t,x)
dx = zeros(3,1);
p7 = t;   %p7=a*b^y*c^z;
p1 = 1; p2 = 2; p3 = 3; p4 = 4; p5 = 5; p6 = 6; p8 = 8;
dx(1) = ((p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))).*x(1));              
dx(2) = ((-((1/p4).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1));
dx(3) = (((-((1/p6).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1))+(p7.*(p8-x(3))));

and call the solver in the command line

[t,x] = ode45('fcn',[0 1],[0 0 0]');

You can set p1 to p8 in the fcn.m file.

[Nivedhitha S]

Sorry i guess you mistook my question. My derivatibve functions are already within a function file:

%solving differential equations

function[dx]=try_eqns(t,x)
p1= p1; % values of parameters wil be given here:
p2= p2;
p3= p3;
p4= p4;
p5= p5;
p6= p6;
p8= p8;

p7=a*y^b*z^c;
% not a function of t or x, y and z are some input variables, a,b,c are parameter values.
% this parameter changes at every time point accoding to changes in inputs

dx = zeros(3,1);
dx(1) = ((p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))).*x(1));              
dx(2) = ((-((1/p4).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1));
dx(3) = (((-((1/p6).*(p1.*((x(2)/(x(2)+p2)).*(x(3)/(x(3)+p3)))))+p5).*x(1))+(p7.*(p8-x(3))));
end

My question was: since p7 changes at every time point (say every 30 second), i'm unable to apply that here. i tried trying for loop, symbolic solving, etc., the solver is considering only the 1st value of parameter 7 and solving the equations. Please help. I want assistance with using ode for changing parameters