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

2 次查看（过去 30 天）
Rose
Rose 2020-11-8
评论： David Wilson 2020-11-8
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...

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by