|
|
Solving a nonlinear system of equations OR finding all roots of a polynomial?
Posted:
Mar 8, 2012 4:59 AM
|
|
Hi,
I'm trying to solve
(a-L)^2/a*2/(a+b)-p=0
(b-U)^2/b*2/(a+b)-p=0
in Visual Basic (!), under the following constraints:
1. 0<p<1 (it's a probability)
2. a>L>0, b>U>0
A possible example is L=U=1, p = 0.022750132: the only solution which satisfies my constraints is [1.1776225370286649805,1.1776225370286649805]. Initially I thought that the system had to have a simple, closed solution, since it comes out of a simple problem with similar triangles. Little did I know that apparently simple math problems, often pose the hardest of challenges :) Luckily, the great folks at sci.math.symbolic (thanks, Martin, Waldek and all the others, for the invaluable help!!) showed me that the system can be reduced to an equation in b:
2*b^4*(p - 1)*(p - 2) + 2*b^3*(L*p*(p - 2) - U*(p^2 - 8*p + 8))+ b^2*(L^2*p^2 + 8*L*U*p + U^2*(p^2 - 16*p + 24))- 4*b*U^2*(L*p + 2*U*(2 - p)) + 2*U^4*(2 - p) = 0
whose solution I can substitute in the following equation for a:
a = (2*b^3*(p - 1)*(p - 2) + 2*b^2*(L*p*(p - 2) - U*(p^2 - 8*p + 8)) + b*(L^2*p^2 + 8*L*U*p + 4*U^2*(5 - 3*p)) - 4*U^2*(L*p + U*(2 - p)))/(U^2*p*(p - 2))
I tried solving f(b)=0 using Newton's method. Unfortunately, it doesn't always converge to the one root satisfying my constraints, out of the four possible. I guess that depends on the starting point: since b must be positive, I choose as a starting point one of the following three, depending on which one is positive:
1) b0=U (which is the solution in the case p=0)
2) b0=(-Sqrt(p) * Sqrt((4 - 2 * p) * U ^ 2 + (4 - 2 * p) * L * U + p * L ^ 2) - (2 - p) * U - p * L) / (3 * p - 2)
3) b0=(Sqrt(p) * Sqrt((4 - 2 * p) * U ^ 2 + (4 - 2 * p) * L * U + p * L ^ 2) + (p - 2) * U - p * L) / (3 * p - 2)
but that's not enough, and sometimes I don't end up with the good root. I guess I should find all real solutions with some methods for finding polynomial roots. Can you point me to a method which a) is simple to code (I'm not exactly an expert in the field, you know), b) is guaranteed to work (since this must go in an automatic process) and c) doesn't rely on complex arithmetics? VBA doesn't support that. However, I could try to implement a complex number class, if that's strictly necessary.
Of course, if you think that I'm completely off and that there's a far better way to solve the initial problem, feel free to show me :)
Best Regards
Sergio Rossi
|
|
