How to solve simultaneous equations using Newton-Raphson's method?

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)

• Why not try $x=y=0$. – copper.hat Jan 2 '13 at 5:29
• 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 2 '13 at 5:33
• @Martin.kv Do you want an algorithm or an application to solve the equations? – Daryl Jan 2 '13 at 6:30
• An algorithm is the first priority and an application would be an added bonus! – Martin Jan 2 '13 at 6: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 2 '13 at 6:41

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.

• I have no idea how to implement this in code... – Martin Jan 2 '13 at 6: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 2 '13 at 8: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 2 '13 at 8:32
• I plan to use C. I just need the numerical roots. – Martin Jan 2 '13 at 8:37
• 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 2 '13 at 8:40

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)