MATLAB Answers

Solving system of ODEs and Plotting

15 views (last 30 days)
Maria Raheb
Maria Raheb on 20 Feb 2021
Commented: Alan Stevens on 21 Feb 2021
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).

Accepted Answer

Alan Stevens
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.
  3 Comments
Alan Stevens
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

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!