DIfferential equations and Matlab

This is a long one; sorry! I have written an Euler's code which works. Now I just have to use it to solve an equation. Here are the instructions:

Use a time-stepping algorithm to build a model for the virus spread in Texas.

The unit of time is measured in days.

Option 1: δ(t) = { 0.03 if t < 90 ; 0.01 if t ≥ 90} γ(t) = {0.97 if t < 90; 0.99 if t ≥ 90}

c(t) = { 100 − 5t/3 if t < 30; 50 if t ≥ 30} q(t) = { 0.0025 − 0.0025 t/60 if t < 30; 0.00125 if t ≥ 30 }

Goal: Apply your Improved Euler code to approximate a solution to the KermackMcKendrick equations over a 180 day period, using the initial parameter values given and N = 1, 800 time steps. Plot each population group (the components of the outbreak vector) on a single, labeled plot.

My questions:

1) is my F(t,v) written correctly? Was I right to use elementwise each operation?

1) how do I call euler to solve this? In other words: how can I write a script to solve this equation when I can't even define what V1 and V2 are? I think I need a refresher in diff eq.

My code: Euler

function y = ieuler(f,y0,t0,tn,N)

h = (tn-t0)/N; 
t = t0:h:tn;
y = zeros(length(y0), N+1);
y(:,1) = y0;

for k = 1:N
    e = y(:,k) + h*f(t(k), y(:,k)) ;
    y(:,k+1) = y(:,k) + .5*h*(f(t(k), y(:,k))+ f(t(k+1),e));
    
end

my attempt at time stepping and plots.

delta = @(t) .03 .* (t<90)+...
    .01 .* (t>=90);
gamma = @(t) 0.97 .*(t<90)+...
    .99 .*(t>=90);
c = @(t) (100-5*t*(1/3)) .* (t<30)+...
    50.* (t>=30);
q = @(t) (.0025-.0025*t*(1/60)) .*(t<30)+...
    .00125.*(t>=30);
beta = @(t) c(t).*q(t);
P = 29000000;
V0 = [P -100 100 0 0]';

F = @(t,V) [((beta.*V1.*V2)./P); ((beta.*V1.*V2)./P) -((delta+gamma)./nu)*V2;...
   (gamma*V2)./nu; (delta.*V2)./nu];


figure(1);
t = linspace(0,180);
plot(t,delta(t),'r', "linewidth",3)
set(gca,'fontsize',16)

figure(2);
t = linspace(0,180);
plot(t,gamma(t),'b', "linewidth",3)
set(gca,'fontsize',16)

figure(3);
t = linspace(0,60);
plot(t,c(t),'g', "linewidth",3)
set(gca,'fontsize',16)

Answer 1

[Jan]

"V1" is the first component of V. So write it as V(1).

There is no need for elementwise operators here, but they operate on scalars also, so this is no problem.

Your beta is a function. So you have to call it by beta(t), not simply beta:

beta = @(t) c(t).*q(t);
F = @(t,V) [(beta(t) * V(1) * V(2)) / P; ...    

Instead of defining V0, provide the initial value as y0 for calling your Euler function.

[Aaron Millan]

this is wonderful! Thank you. I understand now that the first two parts of V are V1 and V2. And, building on your suggestion, I changed beta to beta(t), gamma to gamma(t), and so on. Thank you!