Solving a non-linear system of equations

I keep trying to solve these equations however the error code keeps saying:

"Error using mupadengine/feval_internal Unable to convert to matrix form because the system does not seem to be linear."

How do I fix my code or equations?

clear all; clc
syms a b c d
% a=X_H2 
% b=X_O 
% c=X_H2O 
% d=X_He 
eq1 = ((a^2 * b) / (c^0.6 * d^0.7)) * 1.5 == 0.007484;
eq2 = a + b + c + d == 1;
eq3 = (((2*a) + (2*c)) / (b + c)) == 2;
eq4 = (((2*a) + (2*c)) / d) == 0.857;
eqs = [eq1, eq2, eq3, eq4];
[A,B] = equationsToMatrix([eq1,eq2,eq3,eq4], [a,b,c,d]);

Error using mupadengine/feval_internal
Unable to convert to matrix form because the system does not seem to be linear.

Error in sym/equationsToMatrix (line 61)
T = eng.feval_internal('symobj::equationsToMatrix',eqns,vars);

X = linsolve(A,B)

Answer 1

[Torsten]

syms a b c d
eq1 = ((a^2 * b) / (c^0.6 * d^0.7)) * 1.5 == 0.007484;
eq2 = a + b + c + d == 1;
eq3 = (((2*a) + (2*c)) / (b + c)) == 2;
eq4 = (((2*a) + (2*c)) / d) == 0.857;
sol_a_bcd = solve([eq2 eq3 eq4],[b c d])

sol_a_bcd = struct with fields:
    b: a
    c: 857/2857 - (3714*a)/2857
    d: 2000/2857 - (2000*a)/2857

eq1 = subs(eq1,[b c d],[sol_a_bcd.b sol_a_bcd.c sol_a_bcd.d])

eq1 = 

sol_a = vpasolve(eq1,a)

sol_a = 

0.10638263438209707221391493640957

sol_bcd = subs(sol_a_bcd,a,sol_a)

sol_bcd = struct with fields:
    b: 0.10638263438209707221391493640957
    c: 0.16167129713156859425884491640702
    d: 0.62556343410423726131332521077384

Answer 2

[Walter Roberson]

You cannot. You have unknowns to a fractional power: the system cannot possibly be linear.

You also have a^2*b which again is not a system that can be linear.

equationsToMatrix() is for creating A and b for use with A\b but your system is nonlinear and \ is not appropriate for such a case.

Your equations have multiple solutions. solve() says there are 150 different solutions. Most of them involve complex values. There might possibly only be one real-valued solution... but with the a^2 in there possibly there are two real-valued solutions with a negative in one of them.

format long g
syms a b c d
Q = @(v) sym(v);
eq1 = ((a^2 * b) / (c^Q(0.6) * d^Q(0.7))) * 1.5 == Q(7484)/Q(10)^6;
eq2 = a + b + c + d == 1;
eq3 = (((2*a) + (2*c)) / (b + c)) == 2;
eq4 = (((2*a) + (2*c)) / d) == Q(857)/Q(10)^3;
eqs = [eq1, eq2, eq3, eq4];

N = 100;
for K = 1 : N
   sols(K) = vpasolve(eqs, rand(4,1));
end
sol_a = [sols.a];
sol_b = [sols.b];
sol_c = [sols.c];
sol_d = [sols.d];
solmat = [sol_a; sol_b; sol_c; sol_d].';
scatter(1:4, real(solmat), 'b');

scatter(1:4, imag(solmat), 'r')

solmat(abs(solmat)<1e-10) = 0;
realsol = all(imag(solmat) == 0,2);
reducedmat = double(solmat(realsol,:));
scatter(1:4, reducedmat)

uniquetol(reducedmat, 'ByRows', true)

ans = 1×4

0.106382634382097 0.106382634382097 0.161671297131569 0.625563434104237

[Walter Roberson]

By the way, if you convert the floating point numbers into rationals, then the exact solutions are

syms z
a = 1 - (2857*root(z^150 - (30000*z^140)/2857 + (420000000*z^130)/8162449 - (3640000000000*z^120)/23320116793 + (21840000000000000*z^110)/66625573677601 - (96096000000000000000*z^100)/190349263996906057 + (320320000000000000000000*z^90)/543827847239160604849 - (823680000000000000000000000*z^80)/1553716159562281848053593 + (1647360000000000000000000000000*z^70)/4438967067869439239889115201 + (5839030643839588707499782875119616000000*z^65)/62072053786021932500895566429428141438378939190704737 - (2562560000000000000000000000000000*z^60)/12682128912902987908363202129257 - (9433005886655232160742783320064000000000*z^55)/62072053786021932500895566429428141438378939190704737 + (3075072000000000000000000000000000000*z^50)/36232842304163836454193668483287249 + (15239104824968064880036806656000000000000*z^45)/186216161358065797502686699288284424315136817572114211 - (2795520000000000000000000000000000000000*z^40)/103517230462996080749631310856751670393 - (24618909248736776865972224000000000000000*z^35)/1675945452222592177524180293594559818836231358149027899 + (1863680000000000000000000000000000000000000*z^30)/295748727432779802701696655117739522312801 - (860160000000000000000000000000000000000000000*z^20)/844954114275451896318747343671381815247672457 + (245760000000000000000000000000000000000000000000*z^10)/2414033904484966067782661160869137846162600209649 - 32768000000000000000000000000000000000000000000000/6896894865113548055655062936603126826486548798967193, z, 1)^10)/2000
b = a
d = root(z^150 - (30000*z^140)/2857 + (420000000*z^130)/8162449 - (3640000000000*z^120)/23320116793 + (21840000000000000*z^110)/66625573677601 - (96096000000000000000*z^100)/190349263996906057 + (320320000000000000000000*z^90)/543827847239160604849 - (823680000000000000000000000*z^80)/1553716159562281848053593 + (1647360000000000000000000000000*z^70)/4438967067869439239889115201 + (5839030643839588707499782875119616000000*z^65)/62072053786021932500895566429428141438378939190704737 - (2562560000000000000000000000000000*z^60)/12682128912902987908363202129257 - (9433005886655232160742783320064000000000*z^55)/62072053786021932500895566429428141438378939190704737 + (3075072000000000000000000000000000000*z^50)/36232842304163836454193668483287249 + (15239104824968064880036806656000000000000*z^45)/186216161358065797502686699288284424315136817572114211 - (2795520000000000000000000000000000000000*z^40)/103517230462996080749631310856751670393 - (24618909248736776865972224000000000000000*z^35)/1675945452222592177524180293594559818836231358149027899 + (1863680000000000000000000000000000000000000*z^30)/295748727432779802701696655117739522312801 - (860160000000000000000000000000000000000000000*z^20)/844954114275451896318747343671381815247672457 + (245760000000000000000000000000000000000000000000*z^10)/2414033904484966067782661160869137846162600209649 - 32768000000000000000000000000000000000000000000000/6896894865113548055655062936603126826486548798967193, z, 1)^10
c = d * 1857/1000 - 1