6

I understand how to find roots of a polynomial equation programmatically using Newton-Raphson method as explained here.

How to find the values of $x$ and $y$ from the simultaneous equation given below, using Newton-Raphson method.

Any body knows of any program that can do this?

I would like to extend this for three variables using three equations etc.

$\sin(3x) + \sin(3y) = 0$ (Eq-1)

$\sin(5x) + \sin(5y) = 0$ (Eq-2)

T. Eskin
  • 8,303
Martin
  • 165
  • 1
    Why not try $x=y=0$. – copper.hat Jan 02 '13 at 05:29
  • 1
    Right, that is one solution. i want to find the different values of x and y programmatically and the above equations is just a sample. i am looking for a general way to solve two simultaneous equations programmatically. – Martin Jan 02 '13 at 05:33
  • @Martin.kv Do you want an algorithm or an application to solve the equations? – Daryl Jan 02 '13 at 06:30
  • An algorithm is the first priority and an application would be an added bonus! – Martin Jan 02 '13 at 06:33
  • The basic algorithm is straightforward, $x_{n+1} = x_n-Df(x_n)^{-1} f(x_n)$, but there are many details, as the basic algorithm may not converge if started far from a solution. – copper.hat Jan 02 '13 at 06:41

2 Answers2

18

Newton's method is, provided an initial guess $x_0$ to $f(x)=0$, you just iterate $x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$. In higher dimensions, there is a straightforward analog. So in your case, define $$f\left(\left[ \begin{array}{c} x \\ y \end{array}\right]\right)=\left[\begin{array}{c} f_1(x,y) \\ f_2(x,y) \end{array}\right]=\left[\begin{array}{c} \sin(3x)+\sin(3y) \\ \sin(5x)+\sin(5y) \end{array}\right]$$

so you throw in a vector of size two and your $f$ returns a vector of size two. The derivative is simply the 2x2 Jacobian matrix here $$J=\left[\begin{array}{cc} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} \end{array}\right].$$

The only thing to be careful about is that now you have vector operations. The $f'(x)$ in the denominator is equivalent to inverting the Jacobian matrix and then you have a matrix vector multiply and then a vector subtraction. So the full equation is $$\left[ \begin{array}{c} x_{n+1} \\ y_{n+1} \end{array}\right]=\left[ \begin{array}{c} x_n \\ y_n \end{array}\right]-\left[\begin{array}{cc} \frac{\partial f_1}{\partial x} & \frac{\partial f_1}{\partial y} \\ \frac{\partial f_2}{\partial x} & \frac{\partial f_2}{\partial y} \end{array}\right]^{-1}_{(x_n,y_n)}*f \left(\left[ \begin{array}{c} x_n \\ y_n \end{array}\right]\right)$$

So the Jacobian is inverted and evaluated at the point $(x_n,y_n)$ and then multiplied by $f$. Note that the matrix is being multiplied on the left. And then you can generalize this to any dimension in exactly the same manner.

Fixed Point
  • 7,909
  • 3
  • 31
  • 48
  • I have no idea how to implement this in code... – Martin Jan 02 '13 at 06:59
  • What language are you trying to write this in? How much programming experience do you have? The math here is quite straightforward, especially for the 2x2 case. Are you trying to write a general n-dimensional case or just a 2D case? Are you trying to write for any general function $f$ or just this specific $f$? If you want to use C, then the example you posted can be modified. The modifications of course depend on how general you want the code to be. – Fixed Point Jan 02 '13 at 08:30
  • What is the purpose of this? Do you just need to find a specific solution or are you trying to just learn programming by solving a mathematical problem? Because if you just need the numerical roots, there are other easier ways of doing so. – Fixed Point Jan 02 '13 at 08:32
  • I plan to use C. I just need the numerical roots. – Martin Jan 02 '13 at 08:37
  • 2
    You didn't answer the rest of my questions. Which roots do you need? Because you do have to provide an initial guess and then the method will give you a root...sometimes...and it may not be the closest one to your initial guess...and it may not be the one you want. In order to be helpful, we need a lot more detail from you. And there are infinite roots by the way. Any integral multiple of $\pi$ will work like $x=y=\pi,2\pi,3\pi,4\pi,...$ and so on. – Fixed Point Jan 02 '13 at 08:40
  • Sorry for the delay in replying. The purpose is to do Selective Harmonic Elimination. if the roots are x1, x2, x3 etc the values of x1, x2 need to be between 0 and 90 where x1 < x2 < x3 etc. The solutions needs to be genaralized where there are n unknowns and n equations. – Martin Jan 04 '13 at 05:12
  • https://en.wikibooks.org/wiki/Fractals/Mathematics/Newton_method#systems_of_nonlinear_equations – Adam Oct 07 '16 at 07:35
  • When does this method fails.? – Luai Ghunim Jun 18 '18 at 13:13
  • @LuaiGhunim There are lots of ways Newton's method can fail. It could be that your function is non-smooth and doesn't have a well-defined Jacobian. It could be that the Jacobian is just difficult to approximate because sometimes you don't have the analytic expressions for the partial derivatives and you use finite differences. It could be that the convergence is just really slow. The Jacobian could be near-singular which makes it a pain to invert the Jacobian. The condition number can be enormous and errors can overwhelm your accuracy and precision. – Fixed Point Jun 19 '18 at 16:43
3

The idea behind Newton's method is first-order approximation. If $f(x)$ is differentiable at $x_0$, then near $x_0$ it is similar to the linear function $$ f(x) \sim f(x_0) + (x - x_0) f'(x_0) $$ So if the equation $$ f(x_0) + (x - x_0) f'(x_0) = 0$$ has a solution $x = x_1$, then $f(x_1) \sim 0$. When everything is well-behaved, this will be a better approximation of 0, and repeating the process will converge to a root of $f(x)=0$.

The same idea works in higher dimensions. If we have two functions $f$ and $g$, then

$$f(x,y) \sim f(x_0, y_0) + f_1(x_0, y_0) (x - x_0) + f_2(x_0, y_0) (y - y_0)$$ $$g(x,y) \sim g(x_0, y_0) + g_1(x_0, y_0) (x - x_0) + g_2(x_0, y_0) (y - y_0)$$

If we set the right hand sides to zero and find a solution $(x,y) = (x_1,y_1)$, then $f(x_1,y_1) \sim 0$ and $g(x_1,y_1) \sim 0$.

(Here, $f_1$ means "the derivative of $f$ with respect to its first variable", etc)