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

Question

Rose 2020-11-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/641015-how-to-solve-system-of-odes-with-constant-and-time-varying-coefficients

评论： David Wilson 2020-11-8

MCrain.mat

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)