Solving a non-linear system of equations

16 次查看(过去 30 天)
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)

采纳的回答

Torsten
Torsten 2023-2-6
编辑:Torsten 2023-2-6
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

更多回答(1 个)

Walter Roberson
Walter Roberson 2023-2-6
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
  1 个评论
Walter Roberson
Walter Roberson 2023-2-6
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

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Equation Solving 的更多信息

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by