Initial Parameters are given as which are the same values called upon from (PDE_Coefficients(9)).
Code:
clc; clear; close all;
% Defining Actual Data:clc; clear; close all;
% Defining Actual Data:
global actual
Actual_1D = load('Actual_1D.mat', 'Actual_1D');
Actual_1D = Actual_1D.Actual_1D;
actual = Actual_1D;
% Defining Initial Parameters:
[D1,D2,B_e,B_s,B_a,B_v,theta,lambda,gamma,Omega,rho_e,rho_s,rho_a] = PDE_Coefficients(9);
% Storing Initial Parameters that are being varied:
p0 = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
% lower bound
lb=[10^(-4), 10^(-4), 10^(-4), 10^(-15), 10^(-20), 10^(-8), 10^(-8), 1000, 10^(6), 10^(4)];
% upper bound
ub=[1.0, 1.0, 1.0, 10^(-4), 10^(-4), 10^(-1), 10^(-1), 10^(4), 9*10^(8), 10^(6)];
options = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','Display','iter'); % ,'StepTolerance',10^(-9),'MaxIterations',30 ,'FunctionTolerance',10^(-8)
[p,resnorm,residual,exitflag,output] = lsqnonlin(@PDE_Model,p0,lb,ub,options)
% Plotting:
global t x
global loc
loc = 12; % Specific Location being examined
UU = PDE_Model(p);
figure; hold on; grid on; box on;
plot(t,UU + actual(:,loc))
hold on;
plot(t,actual(:,loc),'ro')
xlabel('t')
ylabel('U')
title(['Population Density at ' num2str(x(loc)) '$$^\circ$$'], 'Interpreter','latex')
legend('Model','Actual Data','Location','best')
figure; hold on;
global Approx
imagesc(x,t,Approx); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Approximated: $$\theta \lambda \hat{E}$$',"Interpreter","latex","FontSize",12)
axis tight
colormap jet
colorbar
hold off;
figure; hold on;
imagesc(x,t,actual); % Contour Plot
grid on; box on;
xlabel('Latitude$^\circ$',"Interpreter","latex","FontSize",12)
ylabel('Time (Weeks)',"Interpreter","latex","FontSize",12)
title('Actual: $$\theta \lambda E$$',"Interpreter","latex",'FontSize',12)
axis tight
colormap jet
colorbar
hold off;
% Printing Table:
% p = [theta, lambda, gamma, B_e, B_v, D1, D2, rho_e, rho_s, rho_a];
fprintf('\n')
fprintf('Optimal Parameter Values \n')
fprintf('-----------------------\n')
fprintf('D1 %3.4e \n', p(6))
fprintf('D2 %3.4e \n', p(7))
fprintf('theta %3.4e \n', p(1))
fprintf('lambda %3.4e \n', p(2))
fprintf('gamma %3.4e \n', p(3))
fprintf('B_e %3.4e \n', p(4))
fprintf('B_v %3.4e \n', p(5))
fprintf('rho_e %3.4e \n', p(8))
fprintf('rho_s %3.4e \n', p(9))
fprintf('rho_a %3.4e \n', p(10))
fprintf('-----------------------\n')
function UU = PDE_Model(p)
% Defining Space and Time:
global A B
[A,B] = Loading_Data(); % Loading Dataframes
global x t t0 tf
% Space (x):
x = sort(A(:,1));
t0 = 9; tf = 37;
% Time (t):
t = [t0:1:tf]';
% Solving Equation:
fprintf('Solving PDE System: \n')
options = odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off','InitialStep',1e-7);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t,options,p);
% Solutions for the 6 Populations:
E = sol(:,:,2);
global actual loc
fprintf('Integrating Newly Infected Symptomatic Data, theta*lambda*E: \n')
N = length(t); theta = p(1); lambda = p(2);
global Approx
Approx = RK4(E,t0,tf,N,theta,lambda);
fprintf('Computing Error Matrix, UU: \n\n')
UU = Approx - actual;
% disp(UU)
UU = UU(:,loc);
% pdepe form:
% c(x,t,u,u_x)u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)
function [c,f,s] = pdefun(x,t,u,dudx,p)
% Parameters:
theta = p(1);
lambda = p(2);
gamma = p(3);
B_e = p(4);
B_v = p(5);
D1 = p(6);
D2 = p(7);
rho_e = p(8);
rho_s = p(9);
rho_a = p(10);
[B_s, B_a] = transmission_constant(t);
Omega = Omega_Surface_Trans(7*t);
% pdepe form:
% c(x,t,u,u_x)*u_t = x^(-m)*(x^(m)*f(x,t,u,u_x))_x + s(x,t,u,u_x)
% c(x,t,u,u_x):
% 6 populations for 6 system of eq'ns
c = ones([6 1]);
% f(x,t,u,u_x,u_y):
f = [D1; D1; D1; D1; D1; D2].*dudx;
% s(x,t,u,u_x):
F1 = -(B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1);
F2 = (B_e*u(2) + B_s*u(3) + B_a*u(4) + B_v*u(6)).*u(1) - lambda*u(2);
F3 = theta*lambda*u(2) - gamma*u(3);
F4 = (1 - theta)*lambda*u(2) - gamma*u(4);
F5 = gamma*(u(3) + u(4));
F6 = rho_e*u(2) + rho_s*u(3) + rho_a*u(4) - Omega*u(6);
F = [F1;F2;F3;F4;F5;F6];
s = F;
end
% Initial Conditions:
function [u0] = pdeic(x,p)
% Input: scalar of space array (x) i.e. 1x1
% Initial Population:
[S0,E0,I_S0,I_A0,R0,V0] = Population(x);
% Storing Initial Pop. into u0:
u0 = [S0; E0; I_S0; I_A0; R0; V0];
end
% Establishing Boundary Conditions:
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t,p)
% Standard Boundary Conditions Solver:
% p(x,t,u) + q(x,t)*f(x,t,u,u_x) = 0
% No Flux Boundary Conditions:
% Domain: (ymin = 32.55, ymax = 33.03)
% Boundary Conditions on the Left:
% S_x(32.55,t) = E_x(32.55,t) = I_s_x(32.55,t) =
% I_a_x(32.55,t) = R_x(32.55,t) = V_x(32.55,t) = 0, etc.
pl = zeros([6 1]);
ql = ones([6 1]);
% Boundary Conditions on the Right:
% S_x(33.03,t) = E_x(33.03,t) = I_s_x(33.03,t) =
% I_a_x(33.03,t) = R_x(33.03,t) = V_x(33.03,t) = 0, etc.
pr = pl;
qr = ql;
end
% Determing B_s and B_a:
function [B_s, B_a] = transmission_constant(t)
% B_s and B_a are defined as such in
% Angleica Bloomquist, Naveen Vaidya, et. al. ”Modeling Risk of Environmental Reservior of SARS-CoV-2”
% Input: t - time (scalar)
% t0 = 9; % initial time
% tf = 38; % final time
tol = 10^(-2); % tolerance condition
n = 1000; % number of points in stage arrays
% First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
t1 = linspace(t0,18,n);
% Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
t2 = linspace(18.01,28,n);
% Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
t3 = linspace(28.01,tf,n);
for i = 1:n
if (t1(i) == t) | (abs(t - t1(i)) <= tol)
% First Stage: 03/07/2020 - May 10, 2020 (~ 64 days ~ 9 weeks)
B_s = 8.92*10^(-8);
B_a = 2.67*10^(-8);
break
elseif (t2(i) == t) | (abs(t - t2(i)) <= tol)
% Second Stage: May 11, 2020 - July 14,2020 (~ 64 days ~ 9 weeks)
B_s = 3.32*10^(-8);
B_a = 2.49*10^(-8);
break
elseif (t3(i) == t) | (abs(t - t3(i)) <= tol)
% Third Stage: July 15, 2020 - September 16, 2020 (~ 63 days ~ 9 weeks)
B_s = 1.64*10^(-8);
B_a = 1.54*10^(-8);
break
else
continue
end
end
end
% Omega function:
function [omega] = Omega_Surface_Trans(t)
Final_tmp = Final_Temperature(t);
% Defining a and b values:
a = 0.094519;
b = -6.4829;
omega = a*Final_tmp + b;
% Note: a and b values are the coefficient averages over all surfaces
end
% Runge-Kutta (4th Order):
function Approx = RK4(E,t0,tf,N,theta,lambda)
h = (tf - t0)/N;
tn = [t0:h:tf]';
% Calculating Cumulative Data:
% f(t) = \theta*\lambda*E(t,x)
% E(t,x) assumed to be 29x30 size
% theta = 0.6;
% lambda = 0.21;
f = @(t,E) theta*lambda*E;
% Initial Condition: y(1) = y0
Y(1,:) = theta*lambda*E(1,:);
% Runge-Kutta Method:
for i = 1:N-1
% k1 = f( t_n, y_n )
k1 = f(tn(i), E(i,:));
% k2 = f( t_n + 0.5*h, y_n + h*0.5*k1 )
k2 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k1);
% k3 = f( t_n + 0.5*h, y_n + h*0.5*k2 )
k3 = f(tn(i) + 0.5*h, E(i,:) + h*0.5*k2);
% k4 = f( t_n + h, y_n + h*k3 )
k4 = f(tn(i) + 0.5*h, E(i,:) + h*k3);
% f(t_n+1, y_n+1)
Y(i+1,:) = Y(i,:) + (1/6)*h*(k1 + 2*k2 + 2*k3 + k4);
end
Approx = Y;
end
function [S0,E0,I_S0,I_A0,R0,V0] = Population(x)
% fprintf('x input: %2.4f \n', x)
[S,E,I_S,I_A,R,V] = Initial_Pop();
tol = 10^(-8); % Defined Tolerance:
% [A,~] = Loading_Data();
Lat = sort(A(:,1)); % Latitude Coordinates:
for i = 1:length(Lat)
if (abs(Lat(i) - x) <= tol)
S0 = S(i);
E0 = E(i);
I_S0 = I_S(i);
I_A0 = I_A(i);
R0 = R(i);
V0 = V(i);
break
else
continue
end
end
end
end
% Calculating Initial Population Function:
function [S0,E0,I_S0,I_A0,R0,V0] = Initial_Pop()
global A B
Lat = A(:,1); Zip = string(A(:,4));
Pop = A(:,3); Case_Data = B(2,:); % Week 1: 03/1/2020 - 03/07/2020
% % Calculating Initial Population Distribution:
% Initial Susceptible Population Density:
[S0] = Initial_Population_Density(Lat,Pop,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
n = length(S0);
% Initial Exposed Population Density:
E0 = zeros(1,n);
% Initial Symptomatic Infected Population Density:
I_S0 = Initial_Population_Density(Lat,Case_Data,Zip,"Initial Susceptible $$S(0,X)$$ Population Density",0);
% Initial Asymptomatic Infected Population Density:
I_A0 = zeros(1,n);
% Initial Recovered Population Density:
R0 = zeros(1,n);
% Initial Environment Population Density:
V0 = zeros(1,n);
end
Final Temperature Function:
% Temperature Functions:
function [F] = Final_Temperature(x)
% Note: x = [0:(1/24):365]' gives the Net Temperature fluctation for One Year
f_s = Seasonal_Temp(x);
f_d = Dirunal_Temp(x);
F = f_s + f_d;
end
% Seasonal Temperature Fluctuation Function:
function [f_s] = Seasonal_Temp(x)
t_0 = 64.712;
epsilon_m = -6.9178;
tau_m = 365;
phi_m = 0.99605;
f_s = t_0 + epsilon_m*sin((2*pi/tau_m)*x + phi_m);
end
% Dirunal Temperature Fluctuation Function:
function [f_d] = Dirunal_Temp(x)
epsilon_d = -4.3372;
tau_d = 1;
phi_d = -5.0872;
m = 0;
f_d = m + epsilon_d*sin((2*pi/tau_d)*x + phi_d);
end
