15 views (last 30 days)

Show older comments

I am trying to graph the solution of a system of ODEs, given by

function system_of_DEs = f(t,x)

p = 0.5; A= 0.7; B= 0.75; C = 0.8;

h = x(1); l = x(2); v = x(3);

dhdt = B*(1-h-l)*h-v*h;

dldt = B*(1-h-l)*l + p*v*h-A*l;

dvdt = A*l+(1-p)*v*h-C*v;

system_of_DEs = [dhdt dldt dvdt]';

end

I want to plot the solution for C<1 for the time range 0-30 and then plot the solution to the same system of DE's on the same plot for C>1 for the time range 30-60.

Is there a way I can do this with one function? (Not have to create two seperate functions for the system of ODEs for C>1 and C<1). I've tried to do that here, but I got an error and I'm not sure why. Here's my function for the C>1 scenario.

function system_of_DEs_one = f(t,x)

p = 0.5; A= 0.7; B= 0.75; C = 1.1;

h = x(1); l = x(2); v = x(3);

dhdt = B*(1-h-l)*h-v*h;

dldt = B*(1-h-l)*l + p*v*h-A*l;

dvdt = A*l+(1-p)*v*h-C*v;

system_of_DEs_one = [dhdt dldt dvdt]';

end

Here's the .m script that I used to solve and plot these systems

timerange= [0 30];

timerange_1 = [30 60];

[t1, hlv1] = ode45(@system_of_DEs, timerange, [0.9,0.1,0.2]);

[t2, hlv2] = ode45(@system_of_DEs_one, timerange_1,hlv1(end,:));

plot(t1, hlv1(:,1), 'LineWidth', 2)

hold on

plot(t2, hlv2(:,1), 'LineWidth', 2)

hold off

I get an error running this and I'm not sure how to debug this.

I would also like if someone could find a way for me to only use one function/define one function for both scenarios described above, like a function that will take C as an input as well- f(t,x,C) and then call it. (I tried doing this and it didn't work--[t1, hlv1] = ode45(@system_of_DEs(1.1), timerange, [0.9,0.1,0.2]); it doesn't work (I put C=1.1, as an example).

Alan Stevens
on 20 Feb 2021

Something like this?

timerange= [0 60];

[t, hvl] = ode45(@f, timerange, [0.9,0.1,0.2]);

plot(t, hvl, 'LineWidth', 2),grid

xlabel('t'),ylabel('h,v,l')

legend('h','v','l')

function system_of_DEs = f(t,x)

p = 0.5; A= 0.7; B= 0.75; C = [0.8 1.1]; %%%% C changed to a vector of values

ptr = 1; if t>30, ptr = 2; end

h = x(1); l = x(2); v = x(3);

dhdt = B*(1-h-l)*h-v*h;

dldt = B*(1-h-l)*l + p*v*h-A*l;

dvdt = A*l+(1-p)*v*h-C(ptr)*v; %%%% C changed to C(ptr)

system_of_DEs = [dhdt dldt dvdt]';

end

Note that one should really use Event detection rather than using the "if" statement within the function. The results look ok to me, but only you can decide if they are acceptable.

Alan Stevens
on 21 Feb 2021

Yes, it is possible to pass extra arguments to the ode function - see below, for example. Try

doc ode45

for more details.

The following also shows one way of plotting the two parts of the curve in different colours.

C = [0.8, 1];

timerange= [0 60];

[t, hvl] = ode45(@(t,x) f(t,x,C), timerange, [0.9,0.1,0.2]);

indx1 = find(t<=30); indx2 = find(t>=30);

t1 = t(indx1); t2 = t(indx2);

h1 = hvl(indx1,1); h2 = hvl(indx2,1); % Similarly for v and l

plot(t1, h1,'r', t2, h2,'b'),grid

xlabel('t'),ylabel('h')

function system_of_DEs = f(t,x,C)

p = 0.5; A= 0.7; B= 0.75;

ptr = 1; if t>30, ptr = 2; end

h = x(1); l = x(2); v = x(3);

dhdt = B*(1-h-l)*h-v*h;

dldt = B*(1-h-l)*l + p*v*h-A*l;

dvdt = A*l+(1-p)*v*h-C(ptr)*v; %%%% C changed to C(ptr)

system_of_DEs = [dhdt dldt dvdt]';

end

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!