Solving system of equations with parameters
이전 댓글 표시
Hi MATLAB community,
I am hoping someone can point me towards some code that I can use to solve a system of 13 equations with 12 unknowns.
All my differential equations are equal to 0 except for the last equation which is the sum of all my states.
State 0 is named S0 etc. Theres a total of 12 states.
I am solving for the values S0, S1, S2, S3, S4 , S5, S6, S7, S8. S9. S10, and S11
Here are the set of 13 equations with 12 unknowns and the model parameters below:
dS0dt = 0 = ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
dS1dt = 0 = ((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
dS2dt = 0 = ((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
dS3dt = 0 = ((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
dS4dt = 0 = ((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
dS5dt = 0 = ((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
dS6dt = 0 = ((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
dS7dt = 0 = ((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
dS8dt = 0 = ((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
dS9dt = 0 = ((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
dS10dt = 0 = ((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
dS11dt = 0 = ((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11);
10000 = S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11;
Ca_sr_conc = 1.3e-3;
MgATP_conc = 5e-3;
MgADP_conc = 36e-6;
Pi_conc = 1e-3;
Ca_cyt_conc = 140e-9;
k_S0_S1 = 4e7;
k_S1_S0 = 450;
k_S1_S2 = 120;
k_S2_S1 = 25;
k_S2_S3 = 1e8;
k_S3_S2 = 16;
k_S3_S4 = 2e7;
k_S4_S3 = 40;
k_S4_S5 = 135;
k_S5_S4 = 765;
k_S5_S6 = 15;
k_S6_S5 = 7.5e4;
k_S6_S7 = 1.5;
k_S7_S6 = 2.0;
k_S7_S8 = 350;
k_S8_S7 = 1.2e5;
k_S8_S9 = 160;
k_S9_S8 = 250;
k_S9_S10 = 300;
k_S10_S9 = 9e4;
k_S10_S11 = 60;
k_S11_S10 = 60;
k_S11_S0 = 250;
k_S0_S11 = 2.7e5;
Thanks in advance for your help!!
댓글 수: 2
Sophie Sophie
2018년 12월 19일
Walter Roberson
2018년 12월 19일
Which MATLAB release are you using? solve() cannot process strings since R2018a.
채택된 답변
madhan ravi
2018년 12월 20일
Like sir Walter says I think you should be using ode solvers to solve the problem because it clearly is a system of odes:
tspan=[0 1]; % time span
ic=[0;1;0;0;0;0;0;0;0;0;0;0]; % initial conditions
[t,x]=ode45(@myod,tspan,ic); % function call
plot(t,x(:,1)) % likewise another 11 solutions can be plotted
function dSdt = myod(t,x) % function definition
S0=x(1);
S1=x(2);
S2=x(3);
S3=x(4);
S4=x(5);
S5=x(6);
S6=x(7);
S7=x(8);
S8=x(9);
S9=x(10);
S10=x(11);
S11=x(12);
Ca_sr_conc = 1.3e-3;
MgATP_conc = 5e-3;
MgADP_conc = 36e-6;
Pi_conc = 1e-3;
Ca=1e-3; % change it to your value
Ca_cyt_conc = 140e-9;
k_S0_S1 = 4e7;
k_S1_S0 = 450;
k_S1_S2 = 120;
k_S2_S1 = 25;
k_S2_S3 = 1e8;
k_S3_S2 = 16;
k_S3_S4 = 2e7;
k_S4_S3 = 40;
k_S4_S5 = 135;
k_S5_S4 = 765;
k_S5_S6 = 15;
k_S6_S5 = 7.5e4;
k_S6_S7 = 1.5; k_S7_S6 = 2.0;
k_S7_S8 = 350;
k_S8_S7 = 1.2e5;
k_S8_S9 = 160;
k_S9_S8 = 250;
k_S9_S10 = 300;
k_S10_S9 = 9e4;
k_S10_S11 = 60;
k_S11_S10 = 60;
k_S11_S0 = 250;
k_S0_S11 = 2.7e5;
dSdt = [ ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11)];
end
추가 답변 (1개)
Walter Roberson
2018년 12월 19일
편집: Walter Roberson
2018년 12월 19일
When you have more equations than unknowns then there is seldom a solution.
It looks like you want Ca to be a model parameter. However, there is no solution of the equations for arbitrary values of Ca, only for one specific value. Ca = -1.84683360434356e-07 -- Yes, it is negative.
syms S0 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11
syms Ca
Ca_sr_conc = 1.3e-3;
MgATP_conc = 5e-3;
MgADP_conc = 36e-6;
Pi_conc = 1e-3;
Ca_cyt_conc = 140e-9;
k_S0_S1 = 4e7;
k_S1_S0 = 450;
k_S1_S2 = 120;
k_S2_S1 = 25;
k_S2_S3 = 1e8;
k_S3_S2 = 16;
k_S3_S4 = 2e7;
k_S4_S3 = 40;
k_S4_S5 = 135;
k_S5_S4 = 765;
k_S5_S6 = 15;
k_S6_S5 = 7.5e4;
k_S6_S7 = 1.5; k_S7_S6 = 2.0;
k_S7_S8 = 350;
k_S8_S7 = 1.2e5;
k_S8_S9 = 160;
k_S9_S8 = 250;
k_S9_S10 = 300;
k_S10_S9 = 9e4;
k_S10_S11 = 60;
k_S11_S10 = 60;
k_S11_S0 = 250;
k_S0_S11 = 2.7e5;
dS0dt = 0 == ((k_S11_S0 * S11) + (k_S1_S0 * S1)) - ((Ca_cyt_conc* k_S0_S1 + Pi_conc * k_S0_S11) * S0);
dS1dt = 0 == ((Ca_cyt_conc* k_S0_S1 * S0) + (k_S2_S1 * S2)) - ((k_S1_S0 + k_S1_S2) * S1);
dS2dt = 0 == ((k_S1_S2 * S1) + (k_S3_S2 * S3)) - ((k_S2_S1 + Ca* k_S2_S3) * S2);
dS3dt = 0 == ((Ca_cyt_conc* k_S2_S3 * S2) + (k_S4_S3 * S4)) - ((k_S3_S2 + MgATP_conc * k_S3_S4) * S3);
dS4dt = 0 == ((MgATP_conc * k_S3_S4 * S3) + (k_S5_S4 * S5)) - ((k_S4_S3 + k_S4_S5) * S4);
dS5dt = 0 == ((k_S4_S5 * S4) + (k_S6_S5 * MgADP_conc * S6)) - ((k_S5_S4 + k_S5_S6) * S5);
dS6dt = 0 == ((k_S5_S6 * S5) + (k_S7_S6 * S7)) - ((k_S6_S5 * MgADP_conc + k_S6_S7) * S6);
dS7dt = 0 == ((k_S6_S7 * S6) + (Ca_sr_conc * k_S8_S7 * S8)) - ((k_S7_S6 * MgADP_conc + k_S7_S8) * S7);
dS8dt = 0 == ((k_S7_S8 * S7) + (k_S9_S8 * S9)) - ((Ca_sr_conc * k_S8_S7 + k_S8_S9) * S8);
dS9dt = 0 == ((k_S8_S9 * S9) + (Ca_sr_conc * k_S10_S9 * S10)) - ((k_S9_S8 + k_S9_S10) * S9);
dS10dt = 0 == ((k_S9_S10 * S9) + (k_S11_S10 * S11)) - ((k_S10_S11 + Ca_sr_conc * k_S10_S9) * S10);
dS11dt = 0 == ((k_S10_S11 * S10) + (Pi_conc * k_S0_S11 * S0)) - ((k_S11_S0 + k_S11_S10) * S11);
tot_SERCA = 10000 == S0 + S1 + S2 + S3 + S4 + S5 + S6 + S7 + S8 + S9 + S10 + S11;
sol = solve([dS0dt,dS1dt,dS2dt,dS3dt,dS4dt,dS5dt,dS6dt,dS7dt,dS8dt,dS9dt,dS10dt,dS11dt,tot_SERCA]);
disp(structfun(@double, sol, 'uniform', 0))
댓글 수: 8
Sophie Sophie
2018년 12월 19일
Walter Roberson
2018년 12월 19일
Your dS2dt involves Ca* k_S2_S3 where Ca has not been defined by you.
You define 12 variables explicitly S0 S1 S2 S3 S4 S5 S6 S7 S8 S9 S10 S11 and you have 13 equations. By declaring Ca to be symbolic that makes it 13 equations in 13 unknowns, and since your equations are linear in the unknowns and your system is not singular, there is a unique solution.
MATLAB permits you to ask to solve for fewer variables than equations, but the answer is almost always empty.
Walter Roberson
2018년 12월 19일
If you leave out the dS2dt and tot_SERCA equation it is possible to solve the remaining set in terms of Ca and S11, the two variables you did not solve for. But then the tot_SERCA reduces to
10000 == (155526453553967425006438293929*S11)*1/4860228527020093611217586850
which forces S11 to a particular value, and the dS2dt reduces to
0 = -375077189670196373016559832188/20250952195917056713406611875*S11 - 235575680000/2349*S11*Ca
and once S11 has been forced, this forces the value for Ca as well.
Sophie Sophie
2018년 12월 19일
편집: Sophie Sophie
2018년 12월 19일
Walter Roberson
2018년 12월 20일
With the revised equation, the solution for the first 12 equations is unique, and it is that all of the variables should be exactly 0. There is no solution in which the sum of them is anything non-zero.
If you have 13 linear equations in 12 variables, the only time there can be a solution is if one or more of the equations are redundant. (Nonlinear equations could potentially have solutions in such a situation though.)
Sophie Sophie
2018년 12월 20일
Sophie Sophie
2018년 12월 20일
Walter Roberson
2018년 12월 20일
The solution of all 0 is correct provided that you have given the correct equations. You can work step by step.
eqn := [0 = 450*S1 - (1378*S0)*1/5 + 250*S11, 0 = (28*S0)*1/5 - 570*S1 + 25*S2,
0 = 120*S1 - 39*S2 + 16*S3, 0 = 14*S2 - 100016*S3 + 40*S4,
0 = 100000*S3 - 175*S4 + 765*S5, 0 = 135*S4 - 780*S5 + (27*S6)*1/10,
0 = 15*S5 - (21*S6)*1/5 + 2*S7,
0 = (3*S6)*1/2 - (6157266382182995*S7)*1/17592186044416 + 156*S8,
0 = 350*S7 - 316*S8 + 250*S9, 0 = 117*S10 - 390*S9,
0 = 300*S9 - 177*S10 + 60*S11, 0 = 270*S0 + 60*S10 - 310*S11];
[ 1378
eqn := [0 = 450 S1 - ---- S0 + 250 S11,
[ 5
28
0 = -- S0 - 570 S1 + 25 S2, 0 = 120 S1 - 39 S2 + 16 S3,
5
0 = 14 S2 - 100016 S3 + 40 S4, 0 = 100000 S3 - 175 S4 + 765 S5,
27 21
0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,
10 5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
]
0 = 300 S9 - 177 S10 + 60 S11, 0 = 270 S0 + 60 S10 - 310 S11]
]
S0_ := solve(eqn[1], S0);
1125 625
S0_ := ---- S1 + --- S11
689 689
eqn2 := simplify(eval(eqn[2 .. -1], S0 = S0_));
[ 386430 3500
eqn2 := [0 = - ------ S1 + ---- S11 + 25 S2,
[ 689 689
0 = 120 S1 - 39 S2 + 16 S3, 0 = 14 S2 - 100016 S3 + 40 S4,
27
0 = 100000 S3 - 175 S4 + 765 S5, 0 = 135 S4 - 780 S5 + -- S6,
10
21
0 = 15 S5 - -- S6 + 2 S7,
5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
303750 44840 ]
0 = ------ S1 - ----- S11 + 60 S10]
689 689 ]
S1_ := solve(eqn2[1], S1);
350 3445
S1_ := ----- S11 + ----- S2
38643 77286
eqn3 := simplify(eval(eqn2[2 .. -1], S1 = S1_));
[ 14000 433459
eqn3 := [0 = ----- S11 - ------ S2 + 16 S3,
[ 12881 12881
0 = 14 S2 - 100016 S3 + 40 S4, 0 = 100000 S3 - 175 S4 + 765 S5,
27 21
0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,
10 5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
786860 253125 ]
0 = - ------ S11 + ------ S2 + 60 S10]
12881 12881 ]
S2_ := solve(eqn3[1], S2);
14000 206096
S2_ := ------ S11 + ------ S3
433459 433459
eqn4 := simplify(eval(eqn3[2 .. -1], S2 = S2_));
[ 196000 43349950000
eqn4 := [0 = ------ S11 - ----------- S3 + 40 S4,
[ 433459 433459
27
0 = 100000 S3 - 175 S4 + 765 S5, 0 = 135 S4 - 780 S5 + -- S6,
10
21
0 = 15 S5 - -- S6 + 2 S7,
5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
26203540 4050000 ]
0 = - -------- S11 + ------- S3 + 60 S10]
433459 433459 ]
S3_ := solve(eqn4[1], S3);
14 433459
S3_ := ------- S11 + ---------- S4
3096425 1083748750
eqn5 := simplify(eval(eqn4[2 .. -1], S3 = S3_));
[ 56000 117048105
eqn5 := [0 = ------ S11 - --------- S4 + 765 S5,
[ 123857 866999
27 21
0 = 135 S4 - 780 S5 + -- S6, 0 = 15 S5 - -- S6 + 2 S7,
10 5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
7487420 3240 ]
0 = - ------- S11 + ------ S4 + 60 S10]
123857 866999 ]
S4_ := solve(eqn5[1], S4);
78400 14738983
S4_ := -------- S11 + -------- S5
23409621 2601069
eqn6 := simplify(eval(eqn5[2 .. -1], S4 = S4_));
[ 392000 13023705 27
eqn6 := [0 = ------ S11 - -------- S5 + -- S6,
[ 867023 867023 10
21
0 = 15 S5 - -- S6 + 2 S7,
5
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
52413380 18360 ]
0 = - -------- S11 + ------ S5 + 60 S10]
867023 867023 ]
S5_ := solve(eqn6[1], S5);
78400 7803207
S5_ := ------- S11 + -------- S6
2604741 43412350
eqn7 := simplify(eval(eqn6[2 .. -1], S5 = S5_));
[ 392000 13056753
eqn7 := [0 = ------ S11 - -------- S6 + 2 S7,
[ 868247 8682470
3 6157266382182995
0 = - S6 - ---------------- S7 + 156 S8,
2 17592186044416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
52486820 16524 ]
0 = - -------- S11 + ------- S6 + 60 S10]
868247 4341235 ]
S6_ := solve(eqn7[1], S6);
3920000 17364940
S6_ := -------- S11 + -------- S7
13056753 13056753
eqn8 := simplify(eval(eqn7[2 .. -1], S6 = S6_));
[ 1960000 26645225141557261584225
eqn8 := [0 = ------- S11 - ----------------------- S7 + 156 S8,
[ 4352251 76565609303995580416
0 = 350 S7 - 316 S8 + 250 S9, 0 = 117 S10 - 390 S9,
0 = 300 S9 - 177 S10 + 60 S11,
263095060 22032 ]
0 = - --------- S11 + ------- S7 + 60 S10]
4352251 4352251 ]
S7_ := solve(eqn8[1], S7);
1379227385882214400 3981411683807770181632
S7_ := ---------------------- S11 + ---------------------- S8
1065809005662290463369 8881741713852420528075
eqn9 := simplify(eval(eqn8[2 .. -1], S7 = S7_));
[ 482729585058775040000
eqn9 := [0 = ---------------------- S11
[ 1065809005662290463369
56525451689785812932020
- ----------------------- S8 + 250 S9, 0 = 117 S10 - 390 S9,
355269668554096821123
7158723507703201131260
0 = 300 S9 - 177 S10 + 60 S11, 0 = - ---------------------- S11
118423222851365607041
6718244744129937408 ]
+ ---------------------- S8 + 60 S10]
2960580571284140176025 ]
S8 := 'S8';
S8_ := solve(eqn9[1], S8);
S8 := S8
24136479252938752000 8881741713852420528075
S8_ := ---------------------- S11 + ---------------------- S9
8478817753467871939803 5652545168978581293202
eqn10 := simplify(eval(eqn9[2 .. -1], S8 = S8_));
[
eqn10 := [0 = 117 S10 - 390 S9, 0 = 300 S9 - 177 S10 + 60 S11, 0 =
[
170849106624052570889180 10077367116194906112
- ------------------------ S11 + ---------------------- S9
2826272584489290646601 2826272584489290646601
]
+ 60 S10]
]
S9_ := solve(eqn10[1], S9);
3
S9_ := -- S10
10
eqn11 := simplify(eval(eqn10[2 .. -1], S9 = S9_));
[
eqn11 := [0 = -87 S10 + 60 S11, 0 =
[
170849106624052570889180 847896891397461486339468 ]
- ------------------------ S11 + ------------------------ S10]
2826272584489290646601 14131362922446453233005 ]
S10_ := solve(eqn11[1], S10);
20
S10_ := -- S11
29
eqn12 := simplify(eval(eqn11[2 .. -1], S10 = S10_));
[ 1563036526507678610428348 ]
eqn12 := [0 = - ------------------------- S11]
[ 81961904950189428751429 ]
S11_ := solve(eqn12[1], S11);
S11_ := 0
And now with S11 being 0, S10 is forced to be 0, and everything else back substitutes to 0. No escaping it.
I would have a question for you, though: the names of your variables suggest that you have differential equations rather than simple linear equations. Are you sure you should be using plain solve() and not dsolve() of something or other?
카테고리
도움말 센터 및 File Exchange에서 Thermal Analysis에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
