Newton-Raphson on a Nonlinear System

Solve two coupled nonlinear equations by iterating the Jacobian-based Newton update, and watch quadratic convergence.

Problem

Newton-Raphson extends the scalar root method to systems: given F(x)=0 with x a vector, iterate x_{k+1} = x_k - J^-1 F(x_k), where J is the Jacobian of first derivatives. Near a simple root convergence is quadratic, meaning the number of correct digits roughly doubles each step.

System

Solve f1 = x^2 + y^2 - 4 = 0 and f2 = x y - 1 = 0. The Jacobian is [[2x, 2y],[y, x]]. A root lies near (1.93, 0.52) where a circle of radius 2 meets the hyperbola xy=1.

import numpy as np
def F(v): x,y=v; return np.array([x*x+y*y-4, x*y-1])
def J(v): x,y=v; return np.array([[2*x,2*y],[y,x]])
v=np.array([2.0,1.0])
for k in range(6):
    dv=np.linalg.solve(J(v),-F(v))
    v=v+dv
    print(k,np.round(v,6),'residual',round(np.linalg.norm(F(v)),2e-1*0+9))

Convergence

Starting from (2,1) the residual norm drops from order 1 to below 1e-9 within four or five iterations, the signature of quadratic convergence. Each step solves a linear system with the Jacobian rather than inverting it explicitly, which is both faster and numerically safer. A poor initial guess can converge to the other root or diverge, so Newton is fast but not globally reliable.

• The Jacobian must be nonsingular at the root; near-singularity slows or breaks convergence.
• When exact derivatives are unavailable, finite-difference or quasi-Newton (Broyden) updates approximate J.
• Kronos equilibrium solvers use Newton iteration on the discretized Grad-Shafranov equation for the plasma boundary.