Need help to solve 5 simultaneous first order differential equations with Initial Condition.

Question

HARSH ZALAVADIYA 2021-4-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/801976-need-help-to-solve-5-simultaneous-first-order-differential-equations-with-initial-condition

评论： HARSH ZALAVADIYA 2021-4-14

I have to solve following first order ordinary differential equations and plot the values of Concentrations (Cmn, Cecm, Crec, Ccirc, Ccells) against time (t).

where value of r(t) is described as following.

Below is my script, but I am getting constant errors. or something i dont understand.

clc; close; clear;
%Constant Parameters
Cmn0 = 6.2E-6; %Initial conc at the MN (
tr = 1805; %release period (s)
ka = 5.01E6; %Association rate (in 1/s)
kd = 5E-4; %Dissociation rate (in 1/s)
ki = 5.05E-3; %Internalization rate (in 1/s)
kc = 5E-3; %Circulation uptake rate (in 1/s)
Rtot = 1.85E-6 ; %initial receptor concentration (in umol/mm^3)
syms r(t) C_mn(t) C_ecm(t) C_rec(t) C_circ(t) C_cells(t) %creating symbolic variable
ode1 = diff(r) == Cmn0/tr;
ode2 = diff(C_mn) == -r;
ode3 = diff(C_ecm) == r- (ka * C_ecm * (Rtot - C_rec - C_cells)) + kd * C_rec - kc * C_ecm;
ode4 = diff(C_rec) == (ka * C_ecm * (Rtot - C_rec - C_cells)) - ((kd + ki)* C_rec);
ode5 = diff(C_circ) == kc * C_ecm;
ode6 = diff(C_cells) == ki * C_rec;
odes = [ode1; ode2; ode3; ode4; ode5; ode6]
cond1 = r(0) == Cmn0/tr;
cond2 = C_mn(0) == 6.2E-6;
cond3 = C_ecm(0) == 0;
cond4 = C_rec(0) == 0;
cond5 = C_circ(0) == 0;
cond6 = C_cells(0) == 0;
conds = [cond1; cond2; cond3; cond4; cond5; cond6];
[VF,Sbs] = odeToVectorField(odes);
odsefcn = matlabFunction(VF,'File', 'Consolvefun')
[t, C] = ode45(@Consolvefun, [0 5000], conds);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2021-4-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/801976-need-help-to-solve-5-simultaneous-first-order-differential-equations-with-initial-condition#answer_675011

在 MATLAB Online 中打开

You have a stiff system (ka~10^6, kd~10^-4), so use ode15s rather than ode45. The following works. I'll leave you to decide if the results make sense!

C0 = [6.2E-6, 0, 0, 0, 0];
tspan = 0:10:5000;
[t, C] = ode15s(@fn, tspan, C0);
C_mn = C(:,1);
C_ecm = C(:,2);
C_rec = C(:,3);
C_circ = C(:,4);
C_cells = C(:,5);
plot(t,C_mn,t,C_ecm,t,C_rec,t,C_circ,t,C_cells),grid
xlabel('t'),ylabel('C')
legend('Cmn','Cecm','Crec','Ccirc','Ccells')
function dCdt = fn(t, C)
    %Constant Parameters
    Cmn0 = 6.2E-6; %Initial conc at the MN (
    tr = 1805; %release period (s)
    ka = 5.01E6; %Association rate (in 1/s)
    kd = 5E-4; %Dissociation rate (in 1/s)
    ki = 5.05E-3; %Internalization rate (in 1/s)
    kc = 5E-3; %Circulation uptake rate (in 1/s)
    Rtot = 1.85E-6 ; %initial receptor concentration (in umol/mm^3)
    r = Cmn0/tr*(t<=tr);
    
    C_mn  = C(1);
    C_ecm = C(2);
    C_rec = C(3);
    C_circ = C(4);
    C_cells = C(5);
    dCdt = [ -r;
             r- ka * C_ecm * (Rtot - C_rec - C_cells) + kd * C_rec - kc * C_ecm;
             ka * C_ecm * (Rtot - C_rec - C_cells) - (kd + ki)* C_rec;
             kc * C_ecm;
             ki * C_rec];
end

3 个评论
显示 1更早的评论隐藏 1更早的评论

HARSH ZALAVADIYA 2021-4-14

And also the condition for r(t) says that when t>tr, it is 0, where as in the graph its not saturating at 1. Thats the part I am confused that how do i put that in code?

HARSH ZALAVADIYA 2021-4-14