0

I am trying to solve two equations simultaneously in Mathematica. The LHS has a high polynomials whereas the RHS has an exponential functions. It does not have to be exact solution, tolerance is needed but I do not how to set them. I am trying to get the value for T1 and T2

This is the equation used to set up function

Th = 1500;
Tc = 300;
F0 = 2.5*10^-4;

\[Epsilon]0 = 0.2;
\[Sigma]0 = 5.670*10^-5; 
R = 2.12*10^-2;
S = 8.82*10^-4;
A0 = 1.2*10^6;
kB = 8.617*10^-5;
wfc = 0.2;
V = 1.0;
Qh = R*\[Sigma]0*(Th^4 - T1^4);
Qc = R*\[Sigma]0*(T2^4 - Tc^4);
I1 = A0*T1^2*Exp[-(wfc + qc*V)/(kB*T1)] ;
I2 = A0*T2^2*Exp[-(wfc)/(kB*T2)];
I = I1 - I2;

D1 = (wfc/qc + V + (2*kB*T1)/qc)*I1 - (wfc/qc + V + (2*kB*T2)/qc)*I2;
D2 = (wfc/qc + (2*kB*T1)/qc)*I1 - (wfc/qc + (2*kB*T2)/qc)*I2;
Dr = \[Epsilon]0*\[Sigma]0 *(T1^4 - T2^4);
Dh = D1 + Dr;
Dc = D2 + Dr;

NSolve[{R*\[Sigma]0 *(Th^4 - T1^4) == F0*Dh, 
  R*\[Sigma]0 *(T2^4 - Tc^4) == F0*Dc, T1 > 0, T2 > 0}, {T1, 
  T2}, Reals]

Below is the msg from mathematica

kamegheka
  • 13
  • 3
  • 1
    Welcome to Mathematica.SE, kamegheka! I suggest the following:
    1. Take the tour and check the faqs.
    2. When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
    3. As you receive help, try to give it too, by answering questions in your area of expertise.
     – Chris K May 26 '19 at 07:54
  • 1
    A few things: 1) I is a built-in symbol, so you can't use it as a variable, hence the Set::wrsym error. It is better to avoid capitalized variables to be safe, unless you know you're using a free symbol. 2) You define Qc but use qc in later definitions. Mathematica is case sensitive so these aren't the same. 3) You will probably need to use FindRoot to solve this equation numerically. Do you have a good initial guess for T1 and T2? – Chris K May 26 '19 at 07:57
  • 1
    If you look at the left and right hand sides of the two equations you will see that you have things like E^11685 and other astonishing quantities that are probably making it very difficult for the solver. You also have the variable qc that doesn't seem to have been given a value. Solving for that too might help. You have I = I1 - I2 when I is defined to be Sqrt[-1] and that will worry readers. Thus far I have not found a way to coax MMA closer to any solution. NMinimize sum of squares cannot find anything close to zero. – Bill May 26 '19 at 08:00

1 Answers1

2

Assume Qc should be renamed qc and then do sum of squares of differences of lhs and rhs

sol=NMinimize[{(R*σ0*(Th^4-T1^4)-F0*Dh)^2+(R*σ0*(T2^4-Tc^4)-F0*Dc)^2, T1> 0&& T2>0}, {T1,T2}]

which promptly returns

{3.6283*^-14, {T1->1504.6217, T2->376.7112}}

Then do this, instead of trying to directly substitute the solution into the original equations because Exp[-111616.] is too small to represent,

{R*σ0 *(Th^4-T1^4), F0*Dh, R*σ0 *(T2^4-Tc^4), F0*Dc}/.sol[[2]]

which instantly returns

{-75346.7562, -75346.7562, 14471.1542, 14471.1542}

You should check the correctness and accuracy of this very carefully. If you change all the approximate decimal constants to exact rationals and use WorkingPrecision->32 in NMinimize then the result is almost the same, but if you push the working precision to 64 then the result is very different, but if you push the working precision to 128 then the result is close to the original again.

Revised to correct my misunderstanding about Qc and qc. My previous method no longer works so use a new method.

Use exact numbers instead of decimal approximations for more precision. Try this without the exact numbers and N[,32] to see that it produces errors and fails. Using other starting points for FindRoot might find other roots.

Th = 1500;Tc = 300;F0 = 25*10^-5;\[Epsilon]0 = 02/10;σ0 = 5670*10^-8;qc=16*10^-20;
R = 212*10^-4;A0 = 12*10^5;kB = 8617*10^-8;wfc = 02/10;V = 1;
I1 = A0*T1^2*Exp[-(wfc + qc*V)/(kB*T1)] ;I2 = A0*T2^2*Exp[-(wfc)/(kB*T2)];
D1 = (wfc/qc + V + (2*kB*T1)/qc)*I1 - (wfc/qc + V + (2*kB*T2)/qc)*I2;
D2 = (wfc/qc + (2*kB*T1)/qc)*I1 - (wfc/qc + (2*kB*T2)/qc)*I2;
Dr = \[Epsilon]0*σ0 *(T1^4 - T2^4);Dh = D1 + Dr;Dc = D2 + Dr;
sol=FindRoot[{R*σ0 *(Th^4 - T1^4)==F0*Dh,R*σ0 *(T2^4 - Tc^4)==F0*Dc},
  {{T1,(300+1500)/2},{T2,(300+1500)/2}},WorkingPrecision->32]
(*{T1->1261.8488582892453396952322322512,T2->1261.8488582892453392649506790059}*)
(*Check result*)
N[{R*σ0 *(Th^4 - T1^4)==F0*Dh,R*σ0 *(T2^4 - Tc^4)==F0*Dc}/.sol,32]
(*{True,True}*)

Please check this very carefully

Bill
  • 12,001
  • 12
  • 13
  • Thanks for your responses. My bad that I did not put qc as a constant. qc (small q) is the charge of an electron which is 1.6x10^-19 coulombs. T1 and T2 should be smaller than 1500 and higher than 300. – kamegheka May 29 '19 at 23:21
  • Hi Bill. I am still having problem in solving the equation. I post a new question on a new problem. Could you please take look at it? Here is the link

    https://mathematica.stackexchange.com/questions/200166/still-having-problem-with-solving-high-polynomials-with-exponential-equations-si

    I try to solve it using NSolve it doesnt have error but doesnt solve it. If I use FindMinimum it says "The functioni value is not a list of real numbers with dimension at (Ta,Tc)={1.,1.}

     – kamegheka Jun 11 '19 at 21:21