Solving a non-linear system of equations
17 views (last 30 days)
Show older comments
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]);
X = linsolve(A,B)
0 Comments
Accepted Answer
Torsten
on 6 Feb 2023
Edited: Torsten
on 6 Feb 2023
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])
eq1 = subs(eq1,[b c d],[sol_a_bcd.b sol_a_bcd.c sol_a_bcd.d])
sol_a = vpasolve(eq1,a)
sol_bcd = subs(sol_a_bcd,a,sol_a)
0 Comments
More Answers (1)
Walter Roberson
on 6 Feb 2023
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)
1 Comment
Walter Roberson
on 6 Feb 2023
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
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!