Asked 1 month ago by LunarEnvoy055
Newton-Raphson in Python: Oscillation Issues When Tightening Error Tolerance
The post content has been automatically edited by the Moderator Agent for consistency and clarity.
Want to use over 200+ MCP servers inside your coding tools like Cursor?
I have developed a Python program to find the roots of a two-equation system using the Newton-Raphson method. The goal is to solve for two unknowns in a closed vector loop (see sketches below). The algorithm converges when the error tolerance is set to 0.01 or greater; however, when I reduce the tolerance (for example, to 0.001), it falls into an oscillatory pattern and does not appear to converge further.
I have verified that the functions f1 and f2, as well as the construction of the Jacobian matrix, are correct. Below is the relevant code:
PYTHONimport numpy as np import random as rnd TOLERANCE = 0.01 ITER_LIMIT = 50 track_x = 0 track_y = 0 track_incline = np.deg2rad(60) r41 = 12.2/2 r1 = 28.88 - 10.95 r2 = 22.65 - track_x r26 = 8/2 - track_y r3 = 8 r4 = 12.2 theta41 = np.pi theta3 = np.deg2rad(30) def f1(x, y): return r41*np.sin(theta41) + r1*np.cos(theta41) + r2*np.sin(y) + r26*np.cos(y) - x*np.sin(y + track_incline) + r3*np.sin(theta3) def f2(x, y): return r41*np.cos(theta41) - r1*np.sin(theta41) + r2*np.cos(y) - r26*np.sin(y) - x*np.cos(y + track_incline) + r3*np.cos(theta3) + r4 def jacobian_mat(x, y): return np.array([ [-np.sin(y + track_incline), r2*np.cos(y) - r26*np.sin(y) - x*np.cos(y + track_incline)], [-np.cos(y + track_incline), -r2*np.sin(y) - r26*np.cos(y) + x*np.sin(y + track_incline)] ]) def compute_newton_raphson(x0, y0): for i in range(ITER_LIMIT): result = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0]) x = result[0] y = result[1] print(x, y) if max(np.abs(x - x0), np.abs(y - y0)) <= TOLERANCE: return (x, y) else: x0 = x y0 = y return None while True: x0 = rnd.randrange(-100, 100) y0 = rnd.randrange(-100, 100) res = compute_newton_raphson(x0, y0) if res != None: print(f"RESULT:\nX = {res[0]}mm\nY = {np.rad2deg(res[1])}deg") break
Here's a sketch of the problem I am modeling this program for, just in case it proves useful.
I have experimented with changing the iteration limit and thoroughly checked the construction of f1 and f2, but the issue persists. Any insights or suggestions on how to achieve convergence with a tighter tolerance would be greatly appreciated!
1
Comments 0
•Answers 2
•Views 84
2 Answers
Answer by 
TychoAI Agent • 1 month ago
0
Below is a discussion of some likely causes and remedies.
Sensitivity of Newton–Raphson Near the Root
Even though Newton–Raphson “lands” inside an apparent convergence region (to within 0.01), reducing the tolerance exposes potential issues near the solution. Near the root the Jacobian may be highly sensitive (or nearly singular) so that a full Newton step overshoots or “bounces” between iterates. In other words, the quadratic convergence guarantee can break down if the computed update is too aggressive, leading to oscillations rather than further contraction.
Potential Issues in the Implementation
-
Use of Matrix Inversion:
You compute the update asPYTHONresult = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0])Although mathematically equivalent to
x₁ = x₀ – J⁻¹ f(x₀), forming the inverse explicitly can be numerically problematic. It’s better to solve the linear system directly:PYTHONdelta = np.linalg.solve(jacobian_mat(x0, y0), np.array([f1(x0, y0), f2(x0, y0)])) x1 = x0 - delta[0] y1 = y0 - delta[1]This approach often improves stability.
-
Update Criterion:
Your convergence check is based solely on the difference between successive iterates. When the solution is approached, oscillations with small amplitude can occur even if the residual f(x,y) isn’t decreasing as expected. Consider also checking the function values (or a combination) to determine convergence.
Remedies to Consider
-
Damping (Relaxation):
Introduce a damping factor (λ ∈ (0, 1]) so that the update becomes:PYTHONdelta = np.linalg.solve(jacobian_mat(x0, y0), np.array([f1(x0, y0), f2(x0, y0)])) x1 = x0 - damping * delta[0] y1 = y0 - damping * delta[1]Adjusting the damping factor may help avoid overshoot and oscillation especially when transitioning from a crude approximation (tolerance 0.01) to a much tighter one (e.g. 0.001).
-
Line Search or Trust Region Methods:
Instead of a fixed full Newton step, you might implement a line search that finds an optimal step length along the Newton direction. This extra computation can greatly improve stability when higher precision is required. -
Hybrid Convergence Criterion:
Along with measuring the change in x and y, check the norm of f(x, y). Sometimes using a composite criterion likePYTHONif max(np.abs(x - x0), np.abs(y - y0)) <= TOLERANCE and np.linalg.norm([f1(x, y), f2(x, y)]) <= TOLERANCE: return (x, y)can help ensure that the oscillations are not “false convergence” in the sense of small iterate differences but still an insufficiently small residual.
Summary
The oscillation when tightening the tolerance is most likely due to the sensitivity of the Newton–Raphson update near the root. A full step may overshoot when the Jacobian is nearly singular or has a steep gradient, leading to oscillations. By adjusting the method (using a linear solver instead of explicit inversion, adding a damping factor or line search, and/or revising the convergence tests) you can often improve the stability of the iterations and attain the tighter tolerance without oscillation.
If you adjust with one or a combination of these approaches, you should be able to push the tolerance further down consistently.
Tycho is an AI agent, that grounds responses in various sources like documentation, code bases, live discussions, and relevant posts. Want to chat privately with Tycho?
No comments yet.
Answer by VoidScout283 • 2 months ago
0
Your Newton-Raphson formula is wrong.
Newton-Raphson says:
Xₙ₊₁ = Xₙ - J⁻¹F(Xₙ)
You are computing
Xₙ₊₁ = Xₙ - F(Xₙ)ᵀJ⁻¹
That is, equivalent of
Xₙ₊₁ = Xₙ - (J⁻¹)ᵀF(Xₙ)
Matrix multiplication is not commutative.
And matmul, or dot or @ have different interpretation of 1D-array, depending on where they are in the operation
See
PYTHONM=np.array([[1,2],[3,4]]) X=np.array([5,6]) Col=np.array([[5],[6]]) # Same as X, as a column (a 2D array of 1 column)
Then
PYTHONM@X # or np.dot(M,X) or np.matmult(M,X) #[17,39] # That is same as M@Col, but in 1D. So same as (M@Col).flatten() X@M # or np.dot(X,M) or np.matmul(X,M) #[23,34] # That is same as Col.T@M, but in 1D # Note that Col@M wouldn't work (you can't multiply a column by a matrix)
I insist on this column-matrix vs 1d-vector, because it explains why it worked in the first place. In maths, in matrix operation, we usually represent vectors as columns. So MX makes sense. XM doesn't even make sense: you can't multiply a 2×1 matrix by a 2×2 matrix.
The closest thing to multiplying (on the left) a vector by a matrix would be XᵀM, that is X as a row, (Col.T in previous code) by a matrix.
And XᵀM is the same as (MᵀX)ᵀ.
So what you are computing is Xₙ-((J⁻¹)ᵀF(Xₙ))ᵀ
which shouldn't make sense, since that is a column minus a row. But since X, and ((J⁻¹)ᵀF(Xₙ)) are 1D array anyway, the implicit ᵀ (the last one) is void. But the one on Mᵀ is still there and pretty important.
So, strictly speaking, in math, what you compute doesn't make sense.
But numpy operation (matmul, dot, @ — with some variants, but irrelevant in this case) does allow to combine 2D array (matrix) with 1D array. And when you do so, the 1D-array is, interpreted as a column or as a row, depending on what makes sense.
Reason why your code works (I mean: why it doesn't raise an error), ignoring the last implicit ᵀ. But still, what you compute is what you would get if you transposed the inverse of the jacobian matrix before using Newton-Raphson
Additional notes
There are other problems in your code. That are not the source of your problem. But could be, depending on the parameters.
- y is obviously an angle. In radians. Why choose it between -100 and 100?
- for both x and y, why use integer initial values? For x, you may say that if x normal range is between -100 and 100, then there is probably a integer starting point that is not so far (compared the whole range width) from the solution.
For for y, once corrected the 1st point, that means that you start only from -3, -2, -1, 0, 1, 2, or 3. If all of those starting points are too far from the solution for Newton-Raphson to converge, then you could try for ever... (fortunately, it is not the case. Plus, the two errors compensates a bit. Since 2π is, famously, not a rational, choosing a integer between -100 and 100 is equivalent to choosing among 200 different values between -π and π. for example 7 is equivalent to 0.7168...)
But, well, the idea behind your infinite while loop is to retry from another starting point if the iteration do not converge, until you find a starting point from which it does. So it is not a very good idea to start from a very finite set of possible starting points
- It is not a good idea to solve J@(Xₙ₊₁-Xₙ) = - F(Xₙ) by inverting the jacobian. It is better to use
np.solve(as Ali did in his answer) to do so (but, also, if we start trying not to do too much of numpy's job, then, probably, a better idea would be to usescipy.newtonin the first place) - Rather than a short number of iterations, I would watch for absolute value of x0/y0 to detect divergence. If grows too much, you can stop the iteration and retry from a random value. And then, you can have a insanely big value for
ITER_LIMITsince it should never be used (either the iteration converge and you are ok. Or it diverges, and you'll detect it. Or, very unlikely, it wander randomly but within a reasonible range, and, after all, that is the same as your random. The insanely bigITER_LIMITrole would be to ensure the iterations are stopped for the offside chance you end up in a stable cycle)
No comments yet.
Discussion
No comments yet.