I am writing a code to solve 5 simultaneous equations with 5 unknowns. I am using the function vpasolve, however the code takes 50 minutes to run. Is there a quicker way of solving the equations?

syms Qa_1 Q1_1 Q2_1 Q3_1 Q4_1;

eqn1 = (Qa_1 == Q1_1 + Q2_1 + Q3_1 + Q4_1);

eqn2 = (Qa_1^2/60.51 + Q1_1^2/0.8616 == 1.035/Q1_1 + 24.3/Qa_1);

eqn3 = (Qa_1^2/60.51 + Q2_1^2/1.346 == 1.321/Q2_1 + 16.57/Qa_1);

eqn4 = (Qa_1^2/60.51 + Q3_1^2/1.346 == 1.236/Q3_1 + 8.873/Qa_1); 

eqn5 = (Qa_1^2/60.51 + Q4_1^2/1.346 == 1.044/Q4_1 + 1.619/Qa_1);

assume (Qa_1, 'real');

assume (Q1_1, 'real');

assume (Q2_1, 'real');

assume (Q3_1, 'real');

assume (Q4_1, 'real');

[sol_Qa_1, sol_Q1_1, sol_Q2_1, sol_Q3_1, sol_Q4_1] = vpasolve([eqn1, eqn2, eqn3, eqn4, eqn5], [Qa_1, Q1_1, Q2_1, Q3_1, Q4_1], [0 Inf; 0 Inf; 0 Inf; 0 Inf; 0 Inf])

Answer 1

[Star Strider]

I would do this numerically, using fsolve. It requires a slight re-write of your equations to make them all implicit.

Example —

syms Qa_1 Q1_1 Q2_1 Q3_1 Q4_1 real
eqn1 = (Qa_1 - (Q1_1 + Q2_1 + Q3_1 + Q4_1));
eqn2 = (Qa_1^2/60.51 + Q1_1^2/0.8616 - (1.035/Q1_1 + 24.3/Qa_1));
eqn3 = (Qa_1^2/60.51 + Q2_1^2/1.346 - (1.321/Q2_1 + 16.57/Qa_1));
eqn4 = (Qa_1^2/60.51 + Q3_1^2/1.346 - (1.236/Q3_1 + 8.873/Qa_1));
eqn5 = (Qa_1^2/60.51 + Q4_1^2/1.346 - (1.044/Q4_1 + 1.619/Qa_1));
Eqnsfcn = matlabFunction([eqn1, eqn2, eqn3, eqn4, eqn5], 'Vars',{[Qa_1, Q1_1, Q2_1, Q3_1, Q4_1]});

B0 = rand(1,5)*100;
[B,fval] = fsolve(Eqnsfcn, B0)

This was almost instantaneous. There are likely multiple roots, so experiment with different initial parameter estimates (here ‘B0’).

EDIT — This version makes it easier to track the individual variable names:

Eqnsfcn = matlabFunction([eqn1, eqn2, eqn3, eqn4, eqn5], 'Vars',{Qa_1, Q1_1, Q2_1, Q3_1, Q4_1});

B0 = rand(1,5)*100;
[B,fval] = fsolve(@(b)Eqnsfcn(b(1),b(2),b(3),b(4),b(5)), B0)

[JS]

Thank you for your help- this speeds up the code a lot. Is there a way of calling one of the results to be used in further equations? For example, using the answer for Qa_1 only in a separate equation later in the code.

[Star Strider]

As always, my pleasure.

To assign the ‘B’ vector elements to their respective variables:

Bc = num2cell(B);
[Qa_1, Q1_1, Q2_1, Q3_1, Q4_1] = Bc{:};

assigning:

Qa_1 =
          6.19400362456342

Q1_1 =
          1.82293308309384

Q2_1 =
          1.91699914162865

Q3_1 =
          1.48226852480147

Q4_1 =
         0.971802875039454

Also, the ‘fval’ output shows that the function very nearly approximated zero for each equation in your system.

[JS]

Thank you, that's very helpful

[Star Strider]

As always, my pleasure.