Clear Filters
Clear Filters

Solve system of simultaneous equations for only real numbers

9 views (last 30 days)
Hello,
I am trying to solve the following system of equations and while I do get an answer, this is complex. For my application, only positive & real solutions are relevant. Furthermore, the solutions to "c1", "c2", "c4" should lie in the range 0 to 1. And c1 + c2 + c3 + c4 should sum to 1.
Is there a way to use vpasolve to find solutions for c1, c2, c4 that are positive, real, and between 0 to 1, or can anyone suggest another solver to use for my problem please?
Thank you.
Code:
clear all;
clc;
a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
syms n c1 c2 c4
init_param = [0; 1];
eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );
eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );
eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );
eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );
sols = vpasolve([eqn1 eqn2 eqn3 eqn4], [n; c1; c2; c4], init_param)
Gives:
Error using sym/vpasolve>checkMultiVarX
Incompatible starting points and variables.
Code:
clear all;
clc;
a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
syms n c1 c2 c4
init_param = [0; 1];
eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );
eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );
eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );
eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );
sols = vpasolve([eqn1 eqn2 eqn3 eqn4], [n; c1; c2; c4], init_param)
Gives:
n: 6.7426863581108857304747603024111 + 3.9194661375046172368785216288436i
c1: 0.030720057908933124309712773265263 + 0.027689129624897200337074934618641i
c2: 0.0022216968838103768366405641539907 + 0.0018149441042006610883472427462662i
c4: 0.017058245207256498853646662580747 - 0.029504073729097861425422177364907i
  2 Comments
Dyuman Joshi
Dyuman Joshi on 8 Mar 2023
Edited: Dyuman Joshi on 8 Mar 2023
Do you know if a real solution exists for all variables?
a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
syms n c1 c2 c4
%4 variablees, need 4 set of initial parameters
init_param = repmat([0 1],4,1)
init_param = 4×2
0 1 0 1 0 1 0 1
eqn1 = ( n == ((log10((c1/a(1,1)) * (a(4,1)/c4))) / (log10(1.57488))) );
eqn2 = ( n == ((log10((c2/a(2,1)) * (a(4,1)/c4))) / (log10(1.55548))) );
eqn3 = ( n == ((log10((c3/a(3,1)) * (a(4,1)/c4))) / (log10(1.30607))) );
eqn4 = ( 1 - (c1 + c2 + c3 + c4) == 0 );
sols = vpasolve([eqn1; eqn2; eqn3; eqn4], [n; c1; c2; c4], init_param)
sols = struct with fields:
n: [0×1 sym] c1: [0×1 sym] c2: [0×1 sym] c4: [0×1 sym]
Grass
Grass on 8 Mar 2023
Hi @Dyuman Joshi, thank you for your answer. I believe that a real solution does exist for all variables... Please see my comment on Fabio's answer. Thank you for the help.

Sign in to comment.

Accepted Answer

Fabio Freschi
Fabio Freschi on 8 Mar 2023
I tried with a numerical solution
clear all;
a = [0.0094; 7.0888e-04; 0.7627; 0.1656];
c3 = 0.95;
% variable mapping
% n -> x(1)
% c1 -> x(2)
% c2 -> x(3)
% c4 -> x(4)
fun = @(x)[log10((x(2)/a(1)) * a(4)/x(4)) / log10(1.57488) - x(1);
log10((x(3)/a(2)) * a(4)/x(4)) / log10(1.55548) - x(1);
log10((c3/a(3)) * a(4)/x(4)) / log10(1.30607) - x(1);
x(2)+x(3)+x(4)+c3-1];
options = optimoptions('fsolve','MaxFunctionEvaluations',10000,'MaxIterations',10000,...
'FunctionTolerance',1e-10);
x = fsolve(fun,[1 1 1 1],options);
The solution has only positive values
>> x
x =
6.9618 0.0431 0.0030 0.0321
however the solution is not very accourate
>> fun(x)
ans =
-0.0006
-0.0000
0.0006
0.0282
Are you sure about the use of parenthesis in your original formulation?
  8 Comments
Alex Sha
Alex Sha on 9 Mar 2023
if taking c3=0.9, there will be one more solution:
n 3.47821500649581
c1 0.021268137459732
c2 0.00153621135446466
c4 0.0771956511857335
if taking c3=0.92, there will be also two solution:
No. 1 2
n 5.49455717047147 8.48765845854466
c1 0.031707727103106 0.055519739381755
c2 0.00223373979164014 0.00376879846933087
c4 0.0460585331052497 0.0207114621489085
Grass
Grass on 9 Mar 2023
Hi @Alex Sha, thanks for this, may I ask how you arrived at these two solutions please? Did you use the same code? Thanks

Sign in to comment.

More Answers (1)

Grass
Grass on 8 Mar 2023
@Dyuman Joshi @Fabio Freschi . I believe you were correct in saying one of the numeric values was incorrect - it seems that if c3 = 0.92 instead of 0.95 - a real, positive soution is found as required. I think I was being too ambitious with what molar composition I could achieve with the distillation. Thanks for the help!

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!