# How to solve this two algebraic equation simultaneously for u and H ?(all real number)

I want to find the non-zero value of H when beta is changing(especially what value of beta when increasing will there be a non-zero H value). H is between [0,100]. This is actually the calculation for a phase-transition temperature.

Thus, I need to find a solution when beta is given a value(like I did in the mathematica code below). But I couldn't even use FindRoot to get a numeric value even for when beta=1.

Need help with solving the two equations below.

A=1;
J=1;
c=1;
beta=1;
FindRoot[{A*beta^(-1.5)*E^(-J*beta-0.5*c*beta*H^2)*(1+Sum[l^(1.5)*E^(-u*l)*Integrate[E^(c*beta*l*H*(1.5x^2-0.5)),{x,0,1}],{l,1,100}])-1==0,
A*beta^(-1.5)*E^(-J*beta-0.5*c*beta*H^2)*(1+Sum[l^(2.5)*E^(-u*l)*Integrate[(1.5x^2-0.5)E^(c*beta*l*H*(1.5x^2-0.5)),{x,0,1}],{l,1,100}])-H==0},{{u,0},{H,10}}]


I have made a serious mistake in my previous question. I missed (1.5x^2-0.5) in the second equation in the integral.

• With the edited formatting, this code outputs {u -> 7.30451, H -> 17.3045}, which suggests to me that you had syntax errors in your original code. Please check your edited post, copy and paste the code, and see if it works. If it does, you can use this to learn proper MMA syntax. – march Aug 31 '15 at 19:47
• @march I get that result as well, but together with a "Failed to converge" error when running the code as is. Do you get no errors? I am on MMA 10.2 Win7-64 – MarcoB Aug 31 '15 at 20:13
• @MarcoB. I do get that error on MMA 10.0 OSX 10.10.5, but for some reason I feel like it worked with no error once. Maybe I'm mis-remembering. Edit: I must be mis-remembering, because that "root" definitely doesn't work when plugged back in. – march Aug 31 '15 at 20:16
• @march I did get that result but when I check it (put the result back to the equation), it clearly is not the solution to the problem. – user56134 Aug 31 '15 at 22:04
• @march When I tried the edited formatting, as you said reaches the result of {u -> 7.30451, H -> 17.3045} But when I put it back to the original equation. I get the value of first equation to be 1.321497032871317*10^368 and value of second equation to be 1.321496435890963*10^370 – user56134 Aug 31 '15 at 22:22

When I first made this answer I was bleary eyed and didn't realize that Willinski's answer matched the question except that he made a nice edit by replacing, for example, 1.5 with 3/2.

This answer is in addition to Willinski's fine work. I followed his procedure.

I wanted to do a numerical study and try to find the region of interest.

A = 1;
J = 1;
c = 1;
beta = 1;

int = Integrate[E^(c*beta*l*H*(1.5 x^2 - 0.5)), {x, 0, 1}]//FullSimplify


sum1 = (1 + Sum[l^(3/2)*E^(-u*l)*int, {l, 1, 100}])

sum2 = (1 + Sum[l^(5/2)*E^(-u*l)*int, {l, 1, 100}])

Note that the integration in sum2 is identical to that used in sum1.

I defined a function representing equation 1, eq1 and equation 2, eq2.

eq1[u_, H_] = A*beta^(-3/2)*E^(-J*beta - 0.5*c*beta*H^2)*(sum1)

eq2[u_, H_] = A*beta^(-3/2)*E^(-J*beta - 0.5*c*beta*H^2)*(sum2)


Note: The function definitions don't use the normal SetDelay (:=) but rather use Set (=).

I first evaluted eq1 and eq2 over different ranges in order to get a birds eye view of where it fit the equation values of 1 and H.

I created data using Table and discriminated the results to lie within the answer for equation 1 (e.g., ≈ 1) and the answer for equation 2 (e.g., ≈ H).

eq1data is for equation 1 and eq2data is for equation 2.

eq1data = Table[{u, H, eq1[u, H]}, {u, 0.1, 2, 0.1}, {H, 0.1, 1, 0.1}];

eq1data =
DeleteCases[
Map[Function[subList, Select[subList, #[[3]] < 2. &]],
eq1data], {}];

eq1data = Partition[Flatten[eq1data], 3];

eq2data = Table[{u, H, eq2[u, H]}, {u, 0.1, 2, 0.1}, {H, 0.1, 1, 0.1}];

eq2data =
DeleteCases[
Map[Function[subList, Select[subList, #[[3]] < 1. &]],
eq2data], {}];

eq2data = Partition[Flatten[eq2data], 3];


I made constant values for each equation (1 for equation 1 and H for equation 2).

eq1Value = Table[{u, H, 1}, {u, 0.7, 2, 0.1}, {H, 0.1, 1, 0.1}];
eq1Value = Partition[Flatten[eq1Value], 3];

eq2Value = Table[{x, y, y}, {x, 0.7, 2, 0.1}, {y, 0.1, 1, 0.1}];
eq2Value = Partition[Flatten[eq2Value], 3];


Finally I plotted the results

Row[{
Show[
ListPlot3D[{eq1Value, eq1data}],
ImageSize -> 400,
AxesLabel -> {u, H, "eq1"}
],
Show[
ListPlot3D[{eq2Value, eq2data}],
ImageSize -> 400,
AxesLabel -> {u, H, "eq2"}
]
}]


The blue area shading is the equation, the orange is the equation value. The first plot is equation 1 and the second plot is equation 2.

You can study the plot and extract values from eq1data and eq2data that show that there is no set of (u,H) values where both equations are satisfied.

If you go very far from the plotted range the equation values grow enormously.

If may well be that with different values for the input parameters there could be a root, but it appears to me that there is no root for the given parameters.

Here is a ContourPlot over the region (very slow to make).

ContourPlot[{eq1[u, H]==1, eq2[u, H]==H}, {u, 0.1, 2}, {H, 0.1, 1},
PlotLegends -> {"eq2", "eq3"}]


This clearly shows that there is no solution for the current parameters in the region of interest.

• Thank you very much for the detailed answer. But I am kind of lost here.data = Table[{u, H, f1[u, H]}, {u, 0.1, 2, 0.1}, {H, 0.1, 1, 0.1}]; data2 = DeleteCases[ Map[Function[subList, Select[subList, #[[3]] < 2. &]], data], {}]; data3 = Partition[Flatten[data2], 3]; Could you give this code more explanation ? – user56134 Sep 1 '15 at 3:24
• How do you know that if {u,H} doesn't have a solution if u,H is outside the region {u, 0.1, 2, 0.1}, {H, 0.1, 1, 0.1} ? – user56134 Sep 1 '15 at 3:29
• I have edited the original answer to follow Willinski's procedure and use nomenclature better suited for the problem. Same result however, with the current parameters there is no solution. In the edited version you no longer have to paste in the sum1* and **sum2 values. You do have to use Set rather than SetDelayed. – Jack LaVigne Sep 1 '15 at 14:44
• Thank you so so much for such an detailed answer to help me. Thank you so much for the help. – user56134 Sep 1 '15 at 17:15
A = 1;
J = 1;
c = 1;
beta = 1;

int = Integrate[E^(c*beta*l*H*(3/2 x^2 - 1/2)), {x, 0, 1}] //FullSimplify;
sum1 = Sum[l^(3/2)*E^(-u*l)*int, {l, 1, 100}];
sum2 = Sum[l^(5/2)*E^(-u*l)*int, {l, 1, 100}];
eq1 = A*beta^(-3/2)*E^(-J*beta - 1/2*c*beta*H^2);
eq2 = eq1*(1 + sum1) - 1;
eq3 = eq1*(1 + sum2) - H;

ContourPlot[{eq2, eq3}, {u, 0, 10}, {H, 0, 5},
PlotLegends -> {"eq2", "eq3"}]


There are no crossings.

You say you have solved for H = 0, but Plot[{eq2, eq3} /. H -> 0, {u, -5, 5}] gives Power::infy: Infinite expression 1/Sqrt[0] encountered. >>

You can yourself increase beta. You don't get crossings in the ContourPlot. For me, the answer to this riddle is found.

• Thank you very much for another detailed answer. I really appreciate this. Helped me so much. – user56134 Sep 1 '15 at 17:26
• One last part, I didn't quite get. How did you be certain that there is no crossing ? There is one part that two lines are really close to each other ? Thank you again. – user56134 Sep 1 '15 at 17:49
• @user56134 Zoom in and you see parallel lines. – user31001 Sep 1 '15 at 18:55
• I have also tried to plot it and zoom in. But I find Sawtooth-like line, and the two line inter-connected with each other @ Willinski – user56134 Sep 1 '15 at 19:03
• @user56134 Please increase the PlotPoints (takes long time for execution). – user31001 Sep 1 '15 at 19:06