How to graph two solutions in one continuous plot?

1 Ansicht (letzte 30 Tage)
Luke Brunot
Luke Brunot am 23 Nov. 2018
Kommentiert: madhan ravi am 24 Nov. 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);

Akzeptierte Antwort

madhan ravi
madhan ravi am 23 Nov. 2018
Bearbeitet: madhan ravi am 23 Nov. 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 Kommentare
Luke Brunot
Luke Brunot am 24 Nov. 2018
Thank you so much madhan ravi! You're the best! That worked!!
madhan ravi
madhan ravi am 24 Nov. 2018
Anytime :) , make sure to accept the answer if it helped

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Mathematics finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by