function Meningitis_Agegroupclear all %This just makes sure the program doesn't remember anything from previous runs!
close all%Declaring the parameter value (Be sure that it follows the same order)
muJ=1/(56*365);pi=muJ*182202000*(1/365);omegaJ=(1/5)*(1/365);kappaJ = 0.0031;nuJ=0.14868;alpha=0.00000986; epsilonJ=((1.90)/2);epsilonA=((1.90)/2); epsilonA1=0;sigmaJ=0.0548;gammaCJ= 0.1118; gammaCJ1= 0.1118*3;gammaIJ=0.1128; gammaIJ1=0.1128*3;deltaJ=0.1923; eta=1; omegaA=(1/5)*(1/365); kappaA=0.0031; nuA=0.14868; nuA1=0;sigmaA=0.0548; gammaCA=0.1118; gammaIA=0.1128;deltaA =0.1923;muA=1/(56*365); betaJ=0.33296/3; betaJ1=0.33296*1.5;betaA=0.33296; betaA1=0.33296*1.5; pars1 = [ muJ pi omegaJ kappaJ nuJ alpha epsilonJ epsilonA sigmaJ gammaCJ gammaIJ deltaJ eta omegaA kappaA nuA sigmaA gammaCA gammaIA deltaA muA betaJ1 betaA1 ]; pars2 = [ muJ pi omegaJ kappaJ nuJ alpha epsilonJ epsilonA sigmaJ gammaCJ1 gammaIJ1 deltaJ eta omegaA kappaA nuA sigmaA gammaCA gammaIA deltaA muA betaJ betaA ]; pars3 = [ muJ pi omegaJ kappaJ nuJ alpha epsilonJ epsilonA1 sigmaJ gammaCJ gammaIJ deltaJ eta omegaA kappaA nuA1 sigmaA gammaCA gammaIA deltaA muA betaJ1 betaA1 ]; pars4 = [ muJ pi omegaJ kappaJ nuJ alpha epsilonJ epsilonA1 sigmaJ gammaCJ1 gammaIJ1 deltaJ eta omegaA kappaA nuA1 sigmaA gammaCA gammaIA deltaA muA betaJ betaA ]; % Set up time vector (times at which you want the solution to the ODE calculated)
t0 = 0;tf = 30;% Define the Initial conditions
SJ0= 5261.332; VJ0=191.908; CJ0=479.77*2; IJ0=152; RJ0=383.816; SA0= 2500; VA0=191.908/2; CA0=479.77/2; IA0=152/2; RA0=383.816/2;x0=[SJ0 VJ0 CJ0 IJ0 RJ0 SA0 VA0 CA0 IA0 RA0];%finally it's time to get a numerical solution!
options = [];[t,y1]=ode45(@Meningitis_Age, [t0,tf], x0, options, pars1); %R0J=3.0213, ROA=0.0048, R0=3.03 | R0>1 Vaccinating both Adults and Juvelines
[t,y2]=ode45(@Meningitis_Age, [t0,tf], x0, options, pars2); %R0J=0.2682, ROA=0.0032, R0=0.27 | R0<1 Vaccinating both Adults and Juvelines
[t,y3]=ode45(@Meningitis_Age, [t0,tf], x0, options, pars3); %R0J=3.0213, ROA=0.0902, R0=3.11 | R0>1 Vaccinating only Juvelines
[t,y4]=ode45(@Meningitis_Age, [t0,tf], x0, options, pars4); %R0J=0.2682, ROA=0.0602, R0=0.33 | R0<1 Vaccinating only Juvelines
set(0,'DefaultAxesFontSize',15) figure(1) plot(t,y1(:,3), 'r', t,y3(:,3), 'g', 'LineWidth',4) legend('Vacc. Both group','Vacc. Only Juveline')xlabel('Time (days)') ylabel('Juveline Carrier Population') figure(2)plot(t,y1(:,8), 'm', t,y3(:,8), 'b', 'LineWidth',4) legend('Vacc. Both group','Vacc. Only Juveline') xlabel('Time (days)') ylabel('Adult Carrier Population') figure(3) plot(t,y1(:,4), 'r', t,y3(:,4), 'g', 'LineWidth',4) legend('Vacc. Both group','Vacc. Only Juveline') xlabel('Time (days)') ylabel('Juveline Infected Population')figure(6)plot(t,y2(:,3), 'r', t,y4(:,3), 'g', 'LineWidth',4)legend('Vacc. Both group','Vacc. Only Juveline')xlabel('Time (days)')ylabel('Juveline Carrier Population')function dydt = Meningitis_Age(t,y,pars)%Declaring my parameters (pars)
muJ = pars(1); pi = pars(2); omegaJ = pars(3); kappaJ = pars(4); nuJ = pars(5); alpha = pars(6); epsilonJ = pars(7); epsilonA = pars(8); sigmaJ = pars(9); gammaCJ = pars(10); gammaIJ = pars(11); deltaJ = pars(12); eta = pars(13); omegaA = pars(14); kappaA = pars(15); nuA = pars(16); sigmaA = pars(17); gammaCA = pars(18); gammaIA = pars(19); deltaA = pars(20); muA =pars(21); betaJ = pars(22); betaA = pars(23); %Declaring the Compartment involved
SJ = max(0,y(1)); VJ= max(0,y(2)); CJ = max(0,y(3)); IJ= max(0,y(4)); RJ = max(0,y(5)); SA = max(0,y(6)); VA= max(0,y(7)); CA = max(0,y(8)); IA= max(0,y(9)); RA = max(0,y(10)); N = SJ + VJ + CJ + IJ + RJ + SA + VA + CA + IA + RA ; lambdaJ = betaJ*(eta*CJ+IJ)/N; lambdaA = betaA*(eta*CA+IA)/N; dSJ = pi+omegaJ*VJ+kappaJ*RJ-(lambdaJ+lambdaA)*SJ-(nuJ+alpha+muJ)*SJ; dVJ = nuJ*SJ-(1-epsilonJ)*(lambdaJ+lambdaA)*VJ-(omegaJ+alpha+muJ)*VJ; dCJ = (lambdaJ+lambdaA)*SJ+(1-epsilonJ)*(lambdaJ+lambdaA)*VJ-(alpha+sigmaJ+gammaCJ+muJ)*CJ; dIJ = sigmaJ*CJ-(alpha+gammaIJ+muJ+deltaJ)*IJ; dRJ = gammaCJ*CJ+gammaIJ*IJ-(alpha+kappaJ+muJ)*RJ; dSA = alpha*SJ+omegaA*VA+kappaA*RA-(lambdaJ+lambdaA)*SA-(nuA+muA)*SA; dVA = alpha*VJ+nuA*SA-(1-epsilonA)*(lambdaJ+lambdaA)*VA-(omegaA+muA)*VA; dCA = alpha*CJ+(lambdaJ+lambdaA)*SA+(1-epsilonA)*(lambdaJ+lambdaA)*VA-(sigmaA+gammaCA+muA)*CA; dIA = alpha*IJ+sigmaA*CA-(gammaIA+muA+deltaA)*IA; dRA = alpha*RJ+gammaCA*CA+gammaIA*IA-(kappaA+muA)*RA; dydt = [dSJ;dVJ;dCJ;dIJ;dRJ;dSA;dVA;dCA;dIA;dRA]; %This is declaring that the function dydt called at the top is the array of [dSJ;dVJ;dCJ;dIJ;dRJ;dSA;dVA;dCA;dIA;dRA] (columnwise)
MATLAB: Hi, how to overcome the error “Vectors must be the same length?”. Am trying to run a simulation, but with maximum time 25 the code works fine, however if I extend the final time it is responding the earlier error stated
dynamical system of odeMATLABwile e. coyote
Related Question
- How to estimate a function to obtain its different values while varying one of its parameter? So that I can obtain the estimate in form of the parameter am varying[nuA], I declared the parameter as a symbol and later declare it as the interval.
- Hello, i need to plot a graph (n) number of times with an offset of a known number between each line on the graph, note all the lnes should be parallel on the same graph but repeated
- Do I get an invalid interpreter on the annotation request
- My annotation of Greek letters fails with capital letters and without them. Any answer
- Exploring Genome-wide Differences in DNA Methylation Profiles documentation
- Shouldn’t there be two parallel lines on the plot
Best Answer