How to graph two solutions in one continuous plot?

조회 수: 1 (최근 30일)
Luke Brunot
Luke Brunot 2018년 11월 23일
댓글: madhan ravi 2018년 11월 24일
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);

채택된 답변

madhan ravi
madhan ravi 2018년 11월 23일
편집: madhan ravi 2018년 11월 23일
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
Luke Brunot
Luke Brunot 2018년 11월 24일
Thank you so much madhan ravi! You're the best! That worked!!
madhan ravi
madhan ravi 2018년 11월 24일
Anytime :) , make sure to accept the answer if it helped

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Mathematics에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by