Changing parameters in an ODE

I wish to change one of the parameters at each time step for an ODE solution plot. My function is as follows:

    function [t,v]=shig(b,p,m,yo)
    [t v] = ode45(@fnsirtry,[0 12],yo);

        function fnsir = fnsirtry(t,v)
            a = 0.25;
            r = 0.14;
            fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
            fnsir(2) = (b*v(1))- (m + r)*v(2);
            fnsir(3) = (r*v(2))-((m+a)*v(3));
            fnsir = fnsir(:);

        end
    end

When I plot the ODE, I want to change the value of b at each step. I have used a code like teh following:

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0000025*temp;
i = [1 2 3 4 5 6 7 8 9 10 11 12];
for i = 1;
    b = beta(:,1)
end
for i = 2;
    b = beta(:,2)
end
for i = 3;
    b = beta(:,3)
end
for i = 4;
    b = beta(:,4)
end
for i = 5;
    b = beta(:,5)
end
for i = 6;
    b = beta(:,6)
end
for i = 7;
    b = beta(:,7)
end
for i = 8;
    b = beta(:,8)
end
for i = 9;
    b = beta(:,9)
end
for i = 10;
    b = beta(:,10)
end
for i = 11;
    b = beta(:,11)
end
for i = 12;
    b = beta(:,12)
end

P = 5000;
m = 0.013;

yo = [200000 160 0]; 
[t v] = shig(b,p,m,yo);
w = [159
148
143
137
102
91
85
137
119
108
104
100];
subplot (1,2,1)
plot(t,v(:,2))
%,'-r*','Linewidth',1.5,'MarkerSize',5)
title('Infected Population')
%legend('Disease free state','Test state')
subplot(1,2,2)
plot(w)

Could I please get some light on this,, pls...

Answer 1

[Star Strider]

I don’t understand all the for loops. Assuming your ODE works and integrates as you want it to (I didn’t run your code), I would just do:

for k1 = 1:12
    [t{k1}, v{k1}] = shig(beta(:,k1),p,m,yo);
end

and then plot the individual cell vectors. Note that if you define your evaluation times as a vector of discrete times for all integrations, rather than as a two-element range, you can use a matrix to store them rather than a cell array. Your call.

[Ojaswita]

Thanks alot for the help. Yes my code runs and integrates appropriately. Could you please explain how to plot individual cell vectors and how to define the evaluation times as a vector.

You have understood what I wanted to do, thanks alot! I tried using your piece of code, but I get the following error:

??? Comma separated list expansion has cell syntax for an array that is not a cell.

[Star Strider]

My pleasure.

First, to use ‘tspan’ (as referred to in the ODE solver documentation) as a vector, just define it as such:

tspan = linspace(0, 12, 50);

creates a 50-element vector going from 0 to 12.

Second, I don’t know precisely what the problem is, but see if changing the loop to this version solves it:

for k1 = 1:12
    [tv, vv] = shig(beta(:,k1),p,m,yo);
    t{k1} = tv;
    v{k1} = vv;
end

The sort of assignment I used initially doesn’t usually cause problems.

Third, to plot individual cell vectors, considering you seem to want to plot them all as subplots, might be something like this:

for k1 = 1:12
    subplot(12, 1, k1)
    plot(t{k1}, v{k1}(:,2))
    grid
end

This is untested code as well, so I leave you to experiment with it to see if it works.

[Ojaswita]

Thanks for the help... the errors still persist on the for loop... I also cant figure out what exactly the problem is as it looks fine to me. I would greatly appreciate if you could try to run my code... I have ti get a separate value of b for each step...

I tried to use parenthesis instead of curly brackets, and I got the following error:

??? In an assignment A(I) = B, the number of elements in B and I must be the same.

[Star Strider]

The parentheses will work if you create a fixed-length time vector:

b = [3 13];                                                 % Proxy For ‘beta’ Array
tv = linspace(0, 2, 75);                                    % Define Fixed-Length Evaluation Times (NECESSARY)
for k1 = 1:2
    ode = @(t,x) [x(1); x(1) - x(2)*sin(b(k1).*pi.*t/2)];   % Define ODE Function With Changing ‘b’-Element
    [t(:,k1),v(:,:,k1)] = ode45(ode, tv, [1; 1]);           % Integrate ODE
end

figure(1)
for k1 = 1:length(b)
    subplot(2,1,k1)
    plot(t(:,k1), v(:,2,k1))
    grid
end

I created a simple ODE for this, but yours should work with it. Be sure to define tspan as a vector and not a range, or this will fail.

I don’t have your constants so I can’t run your code. However this should work with your code, since it does in my simulation. The cell arrays also work in my simulation, so I have no idea why they don’t in your code.

[Ojaswita]

Im getting confused badly... Im trying to put my work in the above format and getting errors... Could you try my code, that would be extremely extremely helpful. My constants are: p = 5000; m =0.013; a =0.25; r = 0.14;

My code (according to the above format) is:

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
tv = linspace(0,12,50);
yo = [200000 160 0]; 
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;

for k1 = 1:12
    ode = @(t,v) [p - m*v(1)-b*v(1)*v(2)+a*v(3); (b*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];    
     [t(:,k1),x(:,:,k1)] = ode45(ode, tv, yo);  
end 

figure(1)
for k1 = 1:length(b)
    subplot(2,1,k1)
    plot(t(:,k1), x(:,2,k1))
    grid
end
-----------------------------------------------------
And when i tried to run the previous assignment format, like this:

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0000025*temp;
P = 5000;
m = 0.013;

yo = [200000 160 0]; 
for k1 = 1:12
    b = beta(:,k1);
    [t{k1}, v{k1}] = shig(b,p,m,yo);
end 

subplot (1,2,1)
plot(t,v(:,2))

I get an error that says:

??? Error using ==> plot Size mismatch in Param Cell / Value Cell pair Error in ==> compare at 64 plot(t,v(:,2))

[Star Strider]

Since you have 12 values of ‘b’, you have to run the integration loop 12 times, each with a different value for ‘b’. Then create 12 subplots.

This works:

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
tv = linspace(0,12,50);
yo = [200000 160 0]; 
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;

for k1 = 1:12
    ode = @(t,v) [p - m*v(1)-b(k1)*v(1)*v(2)+a*v(3); (b(k1)*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];    
     [t(:,k1),x(:,:,k1)] = ode45(ode, tv, yo);  
end 

NrPlots = length(b);
figure(1)
for k1 = 1:NrPlots
    subplot(NrPlots,1,k1)
    plot(t(:,k1), x(:,2,k1))
    grid
end

Note that in this code, I subscripted ‘b’ as b(k1) in your ODE. That, and adjusting the subplots accordingly were the only changes needed.

[Ojaswita]

First of all, I really appreciate your assistance. Thank you so very much!! Its really encouraging!

I do not need subplots... I think we lost track of each other there...

I want to plot the second graph of the solution of my ODE, for twelve points, with each point having a different value of b in the solution.

[Star Strider]

My pleasure.

I’m not quite sure what you mean by ‘twelve points, with each point having a different value of b in the solution.’ The code produces twelve solutions.

See if this does what you want:

figure(2)
plot(t, squeeze(x(:,2,:)))
lgndstr = regexp(sprintf('\\beta = %.3E\n',b), '\n', 'split');
legend(lgndstr(1:end-1), 'Location','SW', 'FontSize',6)
grid

[Ojaswita]

Ok, ill try to be more accurate. I want to change the parameter value of b at every time interval in the ODE solution.

I want to plot the second equation of the ODE solution for 12 time intervals. And at each time interval, I want b to have a different value for the solution.

I hope im more clear...

[Star Strider]

I’m lost. Do you want the ODE solved at 12 points only, changing ‘b’ at each time point? (I thought that’s what I did.)

How do you define the time in that instance? How do you define the times when ‘b’ changes?

If you want to define a time interval and then change ‘b’ at specific times, that’s relatively straightforward. I just need specifically to know what you want.

[Ojaswita]

I tried to further ammend my code as per my understanding but now it appears to be doing the initial thing using only the first value of b.

Here is the function:

function [t,v]=shig(b,p,m,yo)

[t v] = ode45(@fnsirtry,tv,yo);

      function fnsir = fnsirtry(t,v)
          temp = [5 9 12 17 19 24 28 27 22 17 10 7];
          beta = 0.0000025*temp;
          a = 0.25;
          r = 0.14;
          for k=1:12
              b = beta(:,k)
          fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
          fnsir(2) = (b*v(1))- (m + r)*v(2);
          fnsir(3) = (r*v(2))-((m+a)*v(3));
          fnsir = fnsir(:);
         end
      end
  end

And the plotting:

t = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*t
tv = linspace(0,12,50);
yo = [200000 160 0]; 
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;

for k = 1:12
    b(k) = B(:,k)
    [t v]=shig(b(k),p,m,yo)
end 
subplot(1,2,1)
plot(t,v(:,2))

Now it appears to be taking only the first value...

[Star Strider]

Please answer the questions I asked in my previous Comment. I need to understand what you’re doing.

[Ojaswita]

Sorry, I missed one of your comments.

Yes I need the ODE solved at 12 points only, changing ‘b’ at each time point. But I do not want 12 separate subplots, that is what I got. The initial assignment technique made sense to me also, but its giving an error...

b is dependent on temperature, which i have defined as temp. temp gives the average temperature per month. b is a constant value multiplied by the matrix temp.

I need to plot this graph to compare it with original data so that I manage to fix b to give the most perfect match.

[Star Strider]

See if this does what you want:

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
% tv = linspace(0,12,50);
time = 1:12;
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;

for k1 = time
    tv = [k1-0.5, k1, k1+0.5];                  % Three-Element Time Vector
    ode = @(t,v) [p - m*v(1)-b(k1)*v(1)*v(2)+a*v(3); (b(k1)*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];
    [t(:,k1),v(:,:,k1)] = ode45(ode, tv, yo);
    yo = v(3,:,k1);                            % Updated Initial Conditions Vector For Next Iteration
    vp(k1) = yo(2);                            % Element To Be Plotted
end

figure(2)
plot(time(1), vp(1), 'p')
hold on
for k1 = 2:length(time)
    plot(time(k1), vp(k1), 'p')
end
hold off
lgndstr = regexp(sprintf('\\beta = %.3E\n',b), '\n', 'split');
legend(lgndstr(1:end-1), 'Location','SW', 'FontSize',6)
grid

I created a three-element time vector for each integration, and changed ‘b’ with each. The previous end conditions became the new initial conditions. This is the standard way of ‘piecewise-integrating’ an otherwise discontinuous differential equation system. I added ‘vp’ so you could specify it correctly if I got it wrong. It still saves all the results in the ‘v’ matrix.

[Ojaswita]

THANK YOU!!! At least its absolutely somewhere what I needed!

But its still abit different... And im not getting confused in my own mesh... Can we have a conversation out of the forum, please please...

[Star Strider]

My pleasure!

I guess I’m still not understanding exactly what you want to do. What about it is ‘a bit different’ from what you need? How is it different?

I prefer to keep our conversations here because Answers allows markup to make code easier to read, and that may be important. I also can refer to everything else I wrote here, and I know where to look to find all posts in this thread.

[Ojaswita]

I managed to crack my head and get the issue resolved - thank you so much. I was doing a mistake by adding the b definition in both my files hence the code was not working.

But I need some further help with this as the function is complicated. I used a one matrix and an easy formula just to get the code working. Here is the code i used:

function [t,v]=shig(b,p,m,yo)
tv = linspace(0,11,5);
[t v] = ode45(@fnsirtry,tv,yo);
    function fnsir = fnsirtry(t,v)
        a = 0.25;
        r = 0.14;
        fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
        fnsir(2) = (b*v(1))- (m + r)*v(2);
        fnsir(3) = (r*v(2))-((m+a)*v(3));
        fnsir = fnsir(:);
    end
end

and...

temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0025*temp;
tv = linspace(0,12,5);
yo = [200000 160 0]; 
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;

for k = 1:12
    b = beta(k)
    [t,v]=shig(b,p,m,yo)  
end 
subplot(1,2,1)
plot(t,v(:,2))

But I need amendments.. My b is actually a piecewise function.

whereby for 0-6, it is 120*cos(k) for 6-7, it is 50*k - 265 for 7-11 it is 10^k +7

so I need to adjust that.. .would appreciate if you could further help me... And I thank you so much for your help.

[Star Strider]

If it’s piecewise, you have to integrate it over each piece with a different call to ode45.

Do not integrate over the discontinuities! No ode solver does that well.

Use the previous final values of your variables as the initial conditions for the next values of them (or choose a different set, if that works for you). I used the previous values in the code in my comment of 1 Sep 2015 at 17:19, but that is essentially the way to do what you want. It depends on how discontinuous your function is. Note that you have to update your independent variable as well, so it goes from 0-6, 6-7, 7-11. Be sure the end value of the previous interval doesn’t overlap with the first value of the next interval.

[Ojaswita]

Yes, I get it. Thanks so much for the help! I'll do it and come back if I face troubles. Thanks a ton for all the time help! I highly appreciate it. Students like me look forward to such sites for help so that we may learn more and apply better. Thank you!

[Star Strider]

As always, my pleasure!