I WANT TO SOLVE THE FOUR SET OF EQUATIONS TO FIND OUT THE SOLUTION FOR THE FOUR VARIABLES n, p1, p and n2 but couldn't do.

MATLAB CODE

gn1=10^(-8);
gp2=10^(-13);
gn2=10^(-7);
gp1=10^(-6);
E1=0.05;
E2=0.1;
N1=10^15;
N2=10^14;
F=2.5*10^19;
Z=10^15;
T=[100:1:400];
k=8.61*10^(-5);
syms  n p1 p n2
[n,p1,p,n2]= solve(gn1.*n.*p1+ gp1.*p1*F.*exp(-E1./k.*T)== gp1.*p.*( N1-p1), gp2.*p.*n2+ gn2.*n2*F.*exp(-E2./k.*T)== gn2.*n.*( N2-n2), Z== gn1.*n.*p1+ gp2.*p.*n2,((N2-n2)-n)==((N1-p1)-p))

Answer 1

[Walter Roberson]

That is not four equations in four unknowns.

Your T is a vector of length 301, so your first expression in solve() defines a vector of length 301 on the left-hand side, all of which needs to be satisfied by the seeming scalar result of the right-hand side. That has no chance of working out unless all of the 301 left sides have the same result, which could come about if p1 was 0, causing the term involve T to vanish on the left. p1 is in the first term on the left as well, so the left side must be identical to 0, which forces "p" to become 0.

By the same logic, n and n2 must become 0 for the vector in the second expression to work.

With those variables all 0, the right side of the third expression becomes 0. But the left side is non-zero (Z). This is a contradiction, so there is no scalar-valued solution for those symbolic variables.

If you are expecting a vector-valued solution for one of the variables, then create a symbolic vector for it instead of a seeming scalar, and add each of the resulting implied one-to-one equalities in as separate expressions to be satisfied.

There is an alternate approach that might make sense for your purposes: leave T symbolic temporarily, and solve the four equations for [n, n2, p, p1], leaving those in terms of T. You can then evaluate the result over the desired range of T to get vectors of n, n2, p, p1 values. The resulting equation is a polynomial of order 6, and over the range of T from 100 to 400, all of the roots are real-valued. In some cases, some of the variables come out negative, but you did not assert that they needed to be positive (or even real-valued.)

Caution: when you calculate the vectors like above, there is very bad round-off unless you do the calculations at more than 25 digits. 50 digits seems to give answers about as good as 100 digits, but even at 100 digits one of the 6 roots give enormous p2 values suggestive of bad cancellation still being a problem. And calculating even a handful of T values to 50 digits over all 6 roots takes a fair time :(

[Dip Samajdar]

I have tried to solve the problem in the procedure suggested by you. But it's not working.

[Walter Roberson]

Please show code and error message.

In Maple I used (NOTE: symbolic code, not MATLAB code!)

eq1 := [gn1*n*p1+gp1*p1*F*exp(-E1*T/k) = gp1*p*(N1-p1), gp2*p*n2+gn2*n2*F*exp(-E2*T/k) = gn2*n*(N2-n2), Z = gn1*n*p1+gp2*p*n2, N2-n2-n = N1-p1-p];

vars := convert([gn1 = 10^(-8), gp2 = 10^(-13), gn2 = 10^(-7), gp1 = 10^(-6), E1 = 0.5e-1, E2 = .1, N1 = 10^15, N2 = 10^14, F = 2.5*10^19, Z = 10^15, k = 8.61*10^(-5)], rational);

eq1n := eval(eq1, vars);

Q := simplify(solve(eq1n, [n, n2, p, p1]), size):

Then the primary root is

evalf(eval(Q,T=100), 50);

It turns out that to at least a dozen digits of precision, the primary roots will be the same for all T from 100 to 400, because the exp() terms become very small.

For alternate roots, I used Maple's:

allQ := allvalues(Q):
evalf(eval(allQ, T=100), 50);

As mentioned above, this gives 6 roots (including the primary root.) For the first 5 of those roots, the roots will be the same to at least a dozen digits for all T from 100 to 400. However, for the 6th root, the roots will be quite different over T = 100 to 400. For example, for T=100, n2 is -1.0023730696270692585997564216943485205265727045888*10^25229 but for T=400, n2 is -1.1499133814165916393872518188547777254557645418366*10^100890 -- notice the exponent is about 4 times larger. Such values are, of course, not representable in double precision numbers, but a root is a root is a root.

Caution: Q is a pretty long symbolic expression, and allQ is over 6 times larger. Be sure to assign the symbolic form to a variable; if you try to display it to the screen, MuPAD will complain about it being more than 10000 characters. (allQ is about 46500 characters.)