How to graph two solutions in one continuous plot?

1 visualización (últimos 30 días)
Luke Brunot
Luke Brunot el 23 de Nov. de 2018
Comentada: madhan ravi el 24 de Nov. de 2018
I am trying to model the solution to a system of ODEs (written in another file). This file runs and produces a plot, but I can't seem to join the plot for the two different solutions. (There are two ode45 solvers here.) I want the first one to go from t=0 to t=5, the second from t=5 to t=120.
I have tried defining piecewise functions, tried defining two different lin spaces, and it's not working.
Any suggestions would be helpful!
Thank you!
Below is my code:
clear all
clc
axis manual;
close all;
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
lamdax=100;
lamdaz = 600/24;
Dx=.1;
Dy = .5;
Dw = 1/24;
Dz = 1/40;
Dm = 1/1500;
u = 3;
beta = .00061;
eta = .01;
k = 100;
q = 1/30;
c = .006;
Cm = 1.6;
f = 4;
p = .6;
rho = .06;
K1 = 10;
K2 = 10;
K3 = 10;
C1 = .1;
C2 = .1;
C3 = .1;
%parameter values
%solving the ODE @Bliss_ode for time span [0, 1500] and initial
%conditions x1(0)=x2(0)=x3(0)=0.5
[t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
%[t,y] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
%, options);%,[],parfit);
%%plot of x1 versus time
subplot(1,6,1)
plot(t, y(:,1), 'r', 'LineWidth',2);
%%plot of x2 versus time
subplot(1,6,2)
plot(t, y(:,2), 'b', 'LineWidth',2);
%plot of x3 versus time
subplot(1,6,3)
plot(t, y(:,3), 'color',[0.9100 0.4100 0.1700], 'LineWidth',2);
%xlabel('time in days')
%%ylabel('Virus')
subplot(1,6,4)
plot(t, y(:,4), 'k', 'LineWidth',2);
subplot(1,6,5)
plot(t, y(:,5), 'y', 'LineWidth',2);
subplot(1,6,6)
plot(t, y(:,6), 'g', 'LineWidth',2);

Respuesta aceptada

madhan ravi
madhan ravi el 23 de Nov. de 2018
Editada: madhan ravi el 23 de Nov. de 2018
EDITED
clear all
clc
axis manual;
close all;
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
lamdax=100;
lamdaz = 600/24;
Dx=.1;
Dy = .5;
Dw = 1/24;
Dz = 1/40;
Dm = 1/1500;
u = 3;
beta = .00061;
eta = .01;
k = 100;
q = 1/30;
c = .006;
Cm = 1.6;
f = 4;
p = .6;
rho = .06;
K1 = 10;
K2 = 10;
K3 = 10;
C1 = .1;
C2 = .1;
C3 = .1;
%parameter values
%solving the ODE @Bliss_ode for time span [0, 1500] and initial
%conditions x1(0)=x2(0)=x3(0)=0.5
[t,y] = ode45(@measles_ode, [0 150], [lamdax/Dx; .001; 300; lamdaz/Dw; .001; .001]);
[t1,y1] = ode45(@measles_ode, [5, 150], [993.5; 2.75; 300; 595.7; 7.8; .57]);
%, options);%,[],parfit);
%%plot of x1 versus time
subplot(6,1,1)
plot(t, y(:,1), 'r', 'LineWidth',2);
hold on
plot(t1, y1(:,1), 'b', 'LineWidth',2);
%%plot of x2 versus time
subplot(6,1,2)
plot(t, y(:,2), 'b', 'LineWidth',2);
hold on
plot(t1, y1(:,2), 'r', 'LineWidth',2);
%plot of x3 versus time
subplot(6,1,3)
plot(t, y(:,3), 'color',[0.9100 0.4100 0.1700], 'LineWidth',2);
hold on
plot(t1, y1(:,3), 'b', 'LineWidth',2);
%xlabel('time in days')
%%ylabel('Virus')
subplot(6,1,4)
plot(t, y(:,4), 'k', 'LineWidth',2);
hold on
plot(t1, y1(:,4), 'y', 'LineWidth',2);
subplot(6,1,5)
plot(t, y(:,5), 'y', 'LineWidth',2);
hold on
plot(t1, y1(:,5), 'g', 'LineWidth',2);
subplot(6,1,6)
plot(t, y(:,6), 'g', 'LineWidth',2);
hold on
plot(t1, y1(:,6), 'b', 'LineWidth',2);
function dx = measles_ode(t,x)
global lamdax Dx beta q Dy eta k u lamdaz c C1 K1 Dw p C2 K2 rho Dz f Cm C3 K3 Dm
dx = zeros(6,1);
dx(1)= lamdax-Dx*x(1)-beta*q*x(1)*x(3);
dx(2)= beta*q*x(1)*x(3)-Dy*x(2)-eta*x(2)*x(5);
dx(3)= k*x(2)-u*x(3)-beta*q*x(3)*x(1);
dx(4)= lamdaz-((c*q*x(4)*x(3))/(C1*q*x(3)+K1))-Dw*x(4);
dx(5)= ((c*q*x(4)*x(3))/(C1*q*x(3)+K1))+((p*q*x(3)*x(5))/(C2*q*x(3)+K2))-(rho+Dz)*x(5)...
+((f*Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
dx(6)= (rho*x(5))-(Dm*x(6))-((Cm*q*x(3)*x(6))/(C3*q*x(3)+K3));
end
Screen Shot 2018-11-23 at 9.13.18 AM.png
  6 comentarios
Luke Brunot
Luke Brunot el 24 de Nov. de 2018
Thank you so much madhan ravi! You're the best! That worked!!
madhan ravi
madhan ravi el 24 de Nov. de 2018
Anytime :) , make sure to accept the answer if it helped

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Mathematics en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by