# Solve system of simultaneous equations for only real numbers

9 views (last 30 days)
Grass on 8 Mar 2023
Commented: Grass on 9 Mar 2023
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 CommentsShow NoneHide None
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 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.

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?
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 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

### More Answers (1)

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!

### Categories

Find more on Calculus in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!