How to solve system of ODEs with constant and time-varying coefficients?

1 view (last 30 days)
Rose
Rose on 8 Nov 2020
Commented: David Wilson on 8 Nov 2020
I have been trying to solve a system of ODEs in which one coefficient (d1) is a vector and all others are constant. The coefficient d1 is created using some data from MCrain.mat. I have attached the .mat file to this question. I don't know how to implement this situation seamlessly. I have been working on it for many days and have been getting different errors. Right now, the error I am getting is "Sample points must be unique and sorted in ascending order." I looked up online and found that I should use "unique()" but I don't think it would be approapriate since it would change the coefficient d1.
Here is the model equations file where I am constructing d1.
function dxdt = ModelEq_try(t,x,par,tempinit)
dxdt = NaN(7,1);
E = x(1);
L1 = x(2);
L2 = x(3);
L3 = x(4);
L4 = x(5);
P = x(6);
Anext1 = x(7);
% creating a vector of parameters
qe = par(1);
q1 = par(2);
q2 = par(3);
q3 = par(4);
q4 = par(5);
m = par(6);
da = par(7);
%d1=0; If I use d1 to be simply 0, my code runs well. But if I try the following method, i ran into errors.
load('MCrain')
d1 = [];
for time = 1:40*365+1
if MCrain(time)<20
d1(time,1) = 0;
else
d1(time,1) = 0.001;
end
end
d1 = d1';
ft = linspace(1,1,14601);
d1 = interp1(ft,d1,t);
day = t/30;
temp = tempinit + (0.6346*cos(day*0.5163) + 0.7731*sin(day*0.5163));
ne = 0.693;
n1 = 0.693;
n2 = 0.693;
n3 = 0.693;
n4 = 0.693;
np = 0.693;
f1 = 0.666733;
f2 = 0.666733;
f3 = 0.666733;
f4 = 0.666733;
fp = 0.666733-0.056977*temp+ 0.00125224*temp^2;
k = (135000) + 35000*(-1.446*cos(t*0.0172) - 0.7109*sin(t*0.0172)...
-1.347*cos(2*t*0.0172) -1.408*sin(2*t*0.0172) + ...
0.9942*cos(3*t*0.0172) + 0.2396*sin(3*t*0.0172));
dxdt(1) = - ne*E - qe*E;
dxdt(2) = ne*E - (f1 + n1)*L1 - (q1/k)*L1^2 - d1*L1;
dxdt(3) = n1*L1 - (f2 + n2)*L2 - (q2/k)*L2^2;
dxdt(4) = n2*L2 - (f3 + n3)*L3 - (q3/k)*L3^2;
dxdt(5) = n3*L3 - (f4 + n4)*L4 - (q4/k)*L4^2;
dxdt(6) = n4*L4 - (fp + np)*P;
dxdt(7) = np/2*P - (m+da)*Anext1 + m*Anext1;
end
Here is the solve code:
clear all
close all
Q = [0.054,0.054,0.054;...
0.003316,0.000257,0.0008787;...
0.0012,0.0008,0.002;...
0.019,0.0128,0.0088;...
0.0124,0.0229,0.0124;...
];
qe = Q(1,2); q1 = Q(2,2); q2 = Q(3,2); q3 = Q(4,2); q4 = Q(5,2);
%-------------------------------------------------
% paparmeters
%-------------------------------------------------
b = 93.6;
a = 0.049;
households = 4352;
da = 0.167;
di = 0.000255;
dh = 0.000085;
g = 0.48;
m = 1;
epsilon = 1;
% Create vector of parameter values that we will pass to the ode solver:
par = [qe q1 q2 q3 q4 m da];
% Define the initial conditions
E0 = 3954200.8835438923;
L10 = 2310300.5386356956;
L20 = 2007000.7121853685;
L30 = 594100.3681122219;
L40 = 215910.8303970495;
P0 = 103780.6956163843;
A10 = 1289400.7821495022/100/5;
x0 = [E0 L10 L20 L30 L40 P0 A10];
Tf = 40*365;
dt = 1;
t = 1:Tf;
tempinit = 21.25;
sol = ode45(@ModelEq_try,t,x0,[],par,tempinit);
x1 = deval(sol,t);
x = x1';
E = x(:,1);
L1 = x(:,2);
L2 = x(:,3);
plot(t,E)
Any suggestions would be helpful. I am so clueless...
  1 Comment
David Wilson
David Wilson on 8 Nov 2020
Do you really mean the following:
ft = linspace(1,1,14601);
I suspect you mean something like:
ft = linspace(0,1,14601);

Sign in to comment.

Answers (0)

Tags

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!