Message Boards Message Boards

GROUPS:

Solving system of implicit functions: solution satisfies one equation

Posted 3 months ago
606 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?

equations

Attachments:
6 Replies
Posted 3 months ago

Works fine for me. Try restarting the kernel. Maybe this is the problem Clear[T], not clear[T]. WL is case-sensitive.

Posted 3 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!

Posted 3 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.

Posted 3 months ago

That works, many thanks!

Posted 3 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

enter image description here

Posted 3 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]

enter image description here

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]

enter image description here

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract