Solving system of ODEs and Plotting

3 次查看（过去 30 天）
Maria Raheb
Maria Raheb 2021-2-20
评论： Alan Stevens 2021-2-21
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
Alan Stevens 2021-2-20
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 个评论
显示 1更早的评论
Maria Raheb
Maria Raheb 2021-2-20
also I'm guessing that it's not possible to pass arguments into the function the ode45 solves, like ode45(@system_of_DEs(1.1), timerange, [0.9,0.1,0.2]).
Alan Stevens
Alan Stevens 2021-2-21
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by