# Solving system of implicit functions: solution satisfies one equation

Posted 2 months ago
451 Views
|
6 Replies
|
2 Total Likes
|
 I'm trying to solve a system of two equations, one of them is implicit. Using Solve[], I got 2 solutions. However, one of the solutions is not a point for one of the equations, and graphically, the two equations do not cross at the point. I'm confused how Mathematica considers it a solution when it only satisfy one of the equations. Notebook is attached. Here are the two equations, p and T are the only variables.Equation 1: p = \[Phi]/R + ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T)/(R (\[Lambda]/((W - p T) R + \[Phi] T ) + ((1 - \[Lambda]) \[Beta] R)/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))) Equation 2: p = A/T To solve the equations, I did: Solve[p - \[Phi]/R - ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T)/(R (\[Lambda]/((W - p T) R + \[Phi] T ) + ((1 - \[Lambda]) \[Beta] R)/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))) == 0 && p == A/T , {p, T}] and the solutions are: {{p -> -1.98998*10^-8, T -> -5.02518*10^9}, {p -> 0.319279, T -> 313.206}} However, when I substituted T = -5.02518(10^9) into equation 1, I got {p -> 0.325402} instead of {p -> -1.9899810^-8} T -> -5.02518*(10^9) FindRoot[p - \[Phi]/R - ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T)/(R (\[Lambda]/((W - p T) R + \[Phi] T ) + ((1 - \[Lambda]) \[Beta] R)/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))), {p, 0.2}] output is {p -> 0.325402}The parameter values are: \[Beta] = 0.95 R = 1/\[Beta] W = 1000000 A = 100 \[Phi] = 1/5 \[Lambda] = 0.1 Does anyone know how the solution {p -> -0.0000199018, T -> -5.02466*10^6} was derived? It's a point on equation 2 but not equation 1. In the graph below I plotted both equations. p is the y-axis, T is the x-axis. The blue line is p = A/T, why does it look squiggly?  Attachments: Answer
6 Replies
Sort By:
Posted 2 months ago
 Works fine for me. Try restarting the kernel. Maybe this is the problem Clear[T], not clear[T]. WL is case-sensitive. Answer
Posted 2 months ago
 Hello Rohit,Thanks for the tip on Clear[], not clear[]! Now I don't need to restart the kernel each time. My main issue is that one of the solutions ({p -> -1.98998(10^-8), T -> -5.02518(10^9)}) I got from solving the system of equations is not even a point on equation 1. That is, when I plug T = -5.02518(10^9) into equation 1 in my main post, I get {p -> 0.325402}, instead of {p -> -1.9899810^-8} . Would you by any chance know how this can be a solution even when it only satisfies one equation?Thank you! Answer
Posted 2 months ago
 Hi Guangye,Define a symbol for the first equation eq1 = p - \[Phi]/ R + ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T)/(R (\ \[Lambda]/((W - p T) R + \[Phi] T) + ((1 - \[Lambda]) \[Beta] R)/(((W - p T) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))) sol = Solve[eq1 == 0 && p == A/T, {p, T}] (* {{p -> -2.37313*10^-8, T -> -4.21384*10^9}, {p -> 0.0607204, T -> 1646.89}} *) Verify. For eq1 both results are zero to within machine precision. For A / T the result is the two p values from the solution. eq1 /. sol (* {5.55112*10^-17, 0.} *) A / T /. sol (* {-2.37313*10^-8, 0.0607204} *) For an exact solution eq1r = Rationalize@eq1 solr = Solve[eq1r == 0 && p == A/T, {p, T}] Verify as before. Answer
Posted 2 months ago
 That works, many thanks! Answer
Posted 2 months ago
 Hello Rohit,I have a follow-up question. When I plot eq1 around one of the solutions (the negative one), it seems the equation is not actually defined. Would you know why that point would be a solution if it's not on one of the equations? Please see my code and plot below:Thank you again! \[Beta] = 0.95 R = 1/\[Beta] W = 10^6 A = 100 \[Phi] = 0.2 \[Lambda] = 0.1 eq1 = p - \[Phi]/R - ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - T p) R + \[Phi] T) R + \[Phi]*(1 - \[Phi]) T)/(R (\[Lambda]/((W - T p) R + \[Phi] T ) + ((1 - \[Lambda]) \[Beta] R)/(((W - T p) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))) sol = Solve[eq1 == 0 && p == A/T, {p, T}] Which gives 2 solutions: {p -> -0.0000199018, T -> -5.02466(10^6)}, {p -> 0.319279, T -> 313.206}. And then I plot the eq1 for T between -5.025(10^6) and 0 to see what it looks like around the negative solution T= -5.02466(10^6). The plot shows that there is no corresponding value for p when T= -5.02466(10^6). ContourPlot[ p - \[Phi]/R - ((1 - \[Phi]) \[Phi] (1 - \[Lambda]) \[Beta])/(((W - T p) R + \[Phi] T) R + \[Phi]*(1 - \[Phi]) T)/(R (\[Lambda]/((W - T p) R + \[Phi] T ) + ((1 - \[Lambda]) \[Beta] R)/(((W - T p) R + \[Phi] T) R + \[Phi] (1 - \[Phi]) T))) == 0, {T, -5.025*(10^6), 0}, {p, -0.1, 0.4}] T is the x-axis, p is the y-axis  Answer
Posted 2 months ago
 Hi Guangye,By default plots in WL exclude singularities. In this case, they must occur where some denominator in eq1 == 0. Override that. ContourPlot[eq1 == 0, {T, -5.025*(10^6), 0}, {p, -0.1, 0.4}, Exclusions -> None, PlotPoints -> 50] The behavior of the function in that region is quite erratic. Plot3D[eq1, {T, -5.025*10^6, -4*10^6}, {p, -0.1, 0.15}, PlotRange -> All, Exclusions -> None, ColorFunction -> "TemperatureMap", ImageSize -> 600]  Answer