Inverse Kinematics: Working Backward from Where You Want to Be to the Joint Angles

• Introduction — Turning "Put That Cup Right Here" Into Six Angles
• Why Inverse Kinematics Is Harder Than Forward Kinematics
• The Analytic Solution for a 2-Link Arm — Solving It All the Way With the Law of Cosines
• Jacobian — The Matrix That Carries Joint Velocity Into End-Effector Velocity
• Solving Numerically by Iterating With the Pseudoinverse
• Singularities — What Physically Happens
• Damped Least Squares — Surviving Near a Singularity
• Redundancy and Null Space
• Conclusion — Inverse Kinematics Isn't a Problem of Finding an Answer, It's a Problem of Choosing One

Introduction — Turning "Put That Cup Right Here" Into Six Angles

The calculator we built in the forward kinematics post works like this: feed in shoulder 30 degrees, elbow 45 degrees, and it tells you the end effector is at (0.212028, 0.244889). Feed in any angles at all and out comes exactly one answer, done in a handful of multiplications.

But what you actually want to do is the opposite. The cup is at (0.25, 0.15), so put the end effector there — and what angles the joints need to be at is something we have to work out ourselves.

This direction is inverse kinematics (IK), and it's incomparably harder than forward kinematics.

What matters is that the difficulty comes in several different flavors. It's not that the computation is complicated — the nature of the problem itself is different. Some coordinates have no answer at all. Some have two. Some arms have infinitely many. And even where an answer exists, what the neighborhood around that answer looks like differs from spot to spot: in some poses, moving the end effector 1 millimeter only takes 0.15 degrees of joint rotation, while in other poses that same 1 millimeter takes 20.8 degrees.

This post covers those four kinds of difficulty one at a time, and ends with code you can actually use.

Why Inverse Kinematics Is Harder Than Forward Kinematics

Forward kinematics is a function. Feed in a joint-angle vector and exactly one end-effector position comes out. Every point in the domain maps to exactly one value.

Inverse kinematics is the problem of inverting that function — and this function is neither one-to-one nor onto.

Cases with no solution. A point farther than the sum of the link lengths can't be reached by any angle. An arm with a 0.20-meter upper arm and a 0.15-meter forearm has a maximum reach of 0.35 meters. Ask for a point 0.40 meters away and the equation has no real solution. There's also a region on the inside that's unreachable: a point closer than the difference between the two link lengths, 0.05 meters, can't be reached no matter how much you fold the elbow.

Cases with multiple solutions. Most reachable points have exactly two solutions: elbow bent up, or elbow bent down. On a 6-DOF industrial robot, this branching happens at three spots — shoulder, elbow, wrist — producing up to eight solutions.

Cases with infinite solutions. When the joint count exceeds the task's dimensionality, there are continuously infinitely many solutions. Specify only a 2D position with a 3-link planar arm, and the one remaining degree of freedom lets the arm wiggle freely. This is exactly what a human arm does — fix your hand in place and move your elbow up and down, and you can feel it immediately.

Cases with no closed-form solution. Solve a general 6-DOF arm's inverse kinematics algebraically and you get a 16th-degree equation, and a closed-form solution only exists under specific conditions. A structure where the wrist's three axes meet at a single point — a spherical wrist — is the most widely used of those conditions, and it's the practical reason most industrial robots adopt this structure. An arm that doesn't have it has to be solved numerically.

The last row causes the most trouble in practice. When there are two solutions, which one you pick decides the whole arm's pose, and if you solve independently every cycle, the arm can flip instantly the moment the choice switches.

The Analytic Solution for a 2-Link Arm — Solving It All the Way With the Law of Cosines

Let's fully solve the simplest case. Two links, two revolute joints, a 2D position on a plane. Two unknowns, two equations.

Forward kinematics looked like this:

x = L1·cos(θ1) + L2·cos(θ1 + θ2)
y = L1·sin(θ1) + L2·sin(θ1 + θ2)

The key trick is eliminating θ1 first. Square both equations and add them.

x² + y² = L1² + L2² + 2·L1·L2·[cos(θ1)cos(θ1+θ2) + sin(θ1)sin(θ1+θ2)]

The bracketed part is the cosine difference formula, so it becomes cos(θ2).

x² + y² = L1² + L2² + 2·L1·L2·cos(θ2)

This is the law of cosines applied to a triangle. The two links and the line from shoulder to end effector form a triangle, and θ2 is the supplement of that included angle. Rearranged:

D = cos(θ2) = (x² + y² - L1² - L2²) / (2·L1·L2)

Let's plug in numbers. Set the target at (0.212028, 0.244889):

x² + y² = 0.0449559 + 0.0599706 = 0.1049264
L1² + L2² = 0.0400000 + 0.0225000 = 0.0625000
2·L1·L2  = 2 × 0.20 × 0.15 = 0.06

D = (0.1049264 - 0.0625000) / 0.06 = 0.0424264 / 0.06 = 0.7071068

D is 0.7071068. That's cos(45°).

Don't finish here with θ2 = arccos(D). arccos only returns values between 0 and 180 degrees, so it loses the negative solution. Instead, work out the sine as well and feed both into atan2.

sin(θ2) = ±√(1 - D²)
θ2 = atan2(±√(1 - D²), D)

The two signs are exactly the two solutions.

√(1 - 0.7071068²) = √(1 - 0.5) = √0.5 = 0.7071068

θ2 = atan2(+0.7071068, 0.7071068) = +45°
θ2 = atan2(-0.7071068, 0.7071068) = -45°

Now θ1. Take the direction pointing from the shoulder toward the end effector, and subtract the angle the triangle opens up inside it.

θ1 = atan2(y, x) - atan2(L2·sin(θ2), L1 + L2·cos(θ2))

The first term is the angle from the shoulder looking at the end effector, and the second term is how far the upper arm is splayed from that line of sight. Computing with θ2 = +45°:

atan2(0.244889, 0.212028) = 49.1136°
atan2(0.15 × 0.7071068, 0.20 + 0.15 × 0.7071068)
  = atan2(0.1060660, 0.3060660) = 19.1136°

θ1 = 49.1136 - 19.1136 = 30.0000°

Exactly 30 degrees comes out — the same angle we started with in forward kinematics.

Plug in θ2 = -45° and the second term becomes atan2(-0.1060660, 0.3060660) = -19.1136°, so:

θ1 = 49.1136 - (-19.1136) = 68.2271°

The second solution is (68.2271°, -44.9999°). The elbow bends the opposite way, and the end effector lands at the exact same spot.

import numpy as np

L1, L2 = 0.20, 0.15


def ik_2link(x, y, elbow=+1):
    """Analytic solution for a 2-link planar arm. elbow=+1 is elbow-down, -1 is elbow-up."""
    D = (x * x + y * y - L1 * L1 - L2 * L2) / (2 * L1 * L2)
    inner = 1.0 - D * D
    if inner < 0:
        raise ValueError(f"Outside the workspace. D={D:.6f}, 1-D²={inner:.6f}")
    theta2 = np.arctan2(elbow * np.sqrt(inner), D)
    theta1 = np.arctan2(y, x) - np.arctan2(L2 * np.sin(theta2), L1 + L2 * np.cos(theta2))
    return np.array([theta1, theta2])


def fk_2link(q):
    return np.array([L1 * np.cos(q[0]) + L2 * np.cos(q[0] + q[1]),
                     L1 * np.sin(q[0]) + L2 * np.sin(q[0] + q[1])])


target = np.array([0.212028, 0.244889])
for elbow, name in ((+1, "elbow-down"), (-1, "elbow-up  ")):
    q = ik_2link(*target, elbow=elbow)
    print(f"{name}  θ1={np.degrees(q[0]):9.4f}°  θ2={np.degrees(q[1]):9.4f}°  FK check={np.round(fk_2link(q), 6)}")

try:
    ik_2link(0.40, 0.0)
except ValueError as e:
    print("(0.40, 0.00) ->", e)

Here's the output.

elbow-down  θ1=  30.0001°  θ2=  44.9999°  FK check=[0.212028 0.244889]
elbow-up    θ1=  68.2271°  θ2= -44.9999°  FK check=[0.212028 0.244889]
(0.40, 0.00) -> Outside the workspace. D=1.625000, 1-D²=-1.640625

That last line shows exactly how a coordinate outside the workspace shows up in the equations. Plug in (0.40, 0):

D = (0.16 - 0.0625) / 0.06 = 1.625
1 - D² = 1 - 2.640625 = -1.640625

D exceeds 1. Since a cosine can never exceed 1, this means no angle satisfies this target, and the value inside the square root goes negative, halting the computation.

It matters what happens if you skip this check. In numpy, the square root of a negative number isn't an exception — it's nan, and nan propagates silently. θ2 becomes nan, so θ1 becomes nan too, and nan arrives at the function that writes angles to the servo. In C-family languages, converting nan to an integer is undefined behavior, and on an Arduino, if that value reaches servo.write(), the joint jumps to an unpredictable position. The first line of an inverse-kinematics function should always be a reachability check.

Which of the two solutions to pick also has to be decided in code. Three rules are common in practice. One, discard any solution that falls outside the hardware's joint limits. Two, among the remaining solutions, pick the one whose joint-angle change from the current pose is smallest. Three, if there's still more than one, lock in a fixed pose — say, always elbow-up. The second rule matters especially, because without it, the solution can switch mid-trajectory and the arm flips.

Jacobian — The Matrix That Carries Joint Velocity Into End-Effector Velocity

The analytic solution is elegant, but it only works for arms with special structure like a 2-link arm or a spherical wrist. A general arm needs a different approach, and the starting point is the Jacobian.

The idea is this. Differentiate forward kinematics p = f(q):

ṗ = J(q) · q̇

J here is the matrix gathering the partial derivatives of the end-effector coordinates with respect to the joint angles. As components:

J[i][j] = ∂(the i-th component of p) / ∂(the j-th component of q)

A Jacobian is nothing more than a grid of partial derivatives laid out — there isn't some more complicated definition hiding behind it.

Let's derive the 2-link arm's Jacobian by hand. Differentiate x with respect to θ1:

∂x/∂θ1 = -L1·sin(θ1) - L2·sin(θ1 + θ2)
∂x/∂θ2 = -L2·sin(θ1 + θ2)
∂y/∂θ1 =  L1·cos(θ1) + L2·cos(θ1 + θ2)
∂y/∂θ2 =  L2·cos(θ1 + θ2)

Gathered into a matrix:

J = [ -L1·s1 - L2·s12    -L2·s12 ]
    [  L1·c1 + L2·c12     L2·c12 ]

s1 is shorthand for sin(θ1), s12 for sin(θ1+θ2).

Compute this matrix's determinant and a surprisingly simple result comes out.

det(J) = (-L1·s1 - L2·s12)(L2·c12) - (-L2·s12)(L1·c1 + L2·c12)
       = -L1·L2·s1·c12 - L2²·s12·c12 + L1·L2·s12·c1 + L2²·s12·c12
       = L1·L2·(s12·c1 - s1·c12)
       = L1·L2·sin(θ2)

The determinant doesn't depend on θ1 at all — only on θ2. That means how well-behaved the arm is depends entirely on the elbow angle, regardless of where the shoulder points, which makes physical sense — rotating the shoulder just spins the whole arm around.

And when θ2 = 0, the determinant is zero. That's the fully extended pose. The Jacobian becomes singular right there.

Solving Numerically by Iterating With the Pseudoinverse

With the Jacobian in hand, you can turn inverse kinematics into a differential equation and solve it that way.

Say the current joint angles are q, the end effector is at f(q), and the target is p_target. The remaining error is e = p_target - f(q). To eliminate this error, the end effector needs to move by e, and how much the joints need to move to do that is given by J·Δq = e.

If J is square and invertible, Δq = J⁻¹·e. But if the joint count differs from the task dimensionality, it isn't square. That's when you use the Moore-Penrose pseudoinverse J⁺.

When there are more joints than task dimensions (i.e. there's redundancy), J⁺ = Jᵀ(J·Jᵀ)⁻¹, and this solution is the one that minimizes ‖Δq‖. When there are too few joints, J⁺ = (Jᵀ·J)⁻¹Jᵀ, and this solution minimizes the error instead. In both cases, numpy.linalg.pinv handles it automatically.

The whole algorithm looks like this:

1. Pick an initial angle q (usually the current pose)
2. Compute e = p_target - f(q)
3. Stop if ‖e‖ is small enough
4. Δq = J(q)⁺ · e
5. q = q + Δq, go to step 2

Reading it line by line, this is multivariate Newton's method — an iteration searching for the zero of the nonlinear function that is end-effector error.

This method's advantage is that it doesn't care about the arm's structure at all. As long as you can compute the Jacobian, it works identically no matter how many joints there are or in what order they're attached. It has three drawbacks: convergence isn't guaranteed, it lands on different solutions depending on the initial guess, and it blows up near singularities.

That last problem is serious enough to deserve its own section.

Singularities — What Physically Happens

A singularity is a pose where the Jacobian's rank drops. That sounds abstract as a bare definition, but physically, something very concrete happens.

The end effector becomes unable to move in some direction. Extend the arm all the way out and try to push the end effector further away from the shoulder. No matter how you turn any joint, it won't move even a millimeter in that direction — because it's already stretched to the max. This is exactly the situation where the Jacobian's image shrinks from 2D down to 1D.

Singularities on a 6-DOF arm sort into three kinds.

The third one is what you run into most often in practice. When two wrist axes overlap, they end up doing the same job, and the controller starts spinning them furiously in opposite directions to compensate. Most incidents of a wrist suddenly spinning a full turn during teaching on an industrial robot are exactly this.

A singularity isn't binary — it's a matter of continuous degree. The standard measure for that degree is Yoshikawa's manipulability index.

w = √(det(J·Jᵀ))

For a square Jacobian, this equals |det(J)|. An even more practical measure comes from singular value decomposition: the smallest singular value, and the condition number. Let's check this with the 2-link arm we built earlier.

import numpy as np

L1, L2 = 0.20, 0.15


def jacobian(q):
    s1, c1 = np.sin(q[0]), np.cos(q[0])
    s12, c12 = np.sin(q[0] + q[1]), np.cos(q[0] + q[1])
    return np.array([[-L1 * s1 - L2 * s12, -L2 * s12],
                     [ L1 * c1 + L2 * c12,  L2 * c12]])


print(" θ2      det(J)      L1·L2·sin(θ2)    σ1        σ2     condition #")
for d2 in (90, 45, 10, 2, 0):
    J = jacobian(np.radians([30, d2]))
    sigma = np.linalg.svd(J, compute_uv=False)
    cond = sigma[0] / sigma[1] if sigma[1] > 1e-12 else float("inf")
    print(f"{d2:4}°  {np.linalg.det(J):10.7f}  {L1*L2*np.sin(np.radians(d2)):12.7f}"
          f"  {sigma[0]:8.6f}  {sigma[1]:8.6f}  {cond:10.1f}")

Here's the output.

 θ2      det(J)      L1·L2·sin(θ2)    σ1        σ2     condition #
  90°   0.0300000     0.0300000  0.269451  0.111337         2.4
  45°   0.0212132     0.0212132  0.351840  0.060292         5.8
  10°   0.0052094     0.0052094  0.379341  0.013733        27.6
   2°   0.0010470     0.0010470  0.380731  0.002750       138.5
   0°  -0.0000000     0.0000000  0.380789  0.000000         inf

The second and third columns agree on every row. Our hand-derived det(J) = L1·L2·sin(θ2) was correct. The -0.0000000 on the last row is floating-point's signed zero — its value is zero.

What matters is the right three columns. σ1 stays roughly the same, between 0.27 and 0.38, as the elbow angle changes, but σ2 drops from 0.111 all the way to 0. The condition number goes from 2.4 to infinity.

A small singular value is the "gear ratio" in that direction. σ2 = 0.00275 means moving the end effector 1 meter in that direction requires the joints to turn 1/0.00275 = 364 radians. For 1 millimeter, that's 0.364 radians, or 20.8 degrees. In the other direction (σ1 = 0.3807), 1 millimeter only takes 0.0026 radians, or 0.15 degrees. The same 1 millimeter differs by 139x depending on direction.

Damped Least Squares — Surviving Near a Singularity

Now it's clear why the iterative method from earlier blows up. In Δq = J⁺·e, J⁺ contains the reciprocals of the singular values, and as a singular value goes to zero, its reciprocal diverges.

Let's look at the numbers. At θ2 = 1°, command the end effector to move 1 millimeter in the x direction:

|Δq| given by the pseudoinverse = 0.626520 rad = 35.8970°

The joints have to turn 35.9 degrees for that 1 millimeter. At a 100Hz control rate, that's a command to get it done in 10 milliseconds — 3590 degrees per second. A typical hobbyist servo's no-load speed is around 400 degrees per second, so this is nine times over. In practice, the servo pushes at max speed, completely misses the trajectory, and either hits a joint limit or the power supply collapses and the board resets.

The fix is to give up a bit of accuracy. Add a small value to the diagonal when computing the inverse.

Δq = Jᵀ·(J·Jᵀ + λ²·I)⁻¹·e

This is damped least squares (DLS) — known in robotics from Nakamura, Hanafusa, and Wampler's 1986 work, and known in general numerical analysis as the Levenberg-Marquardt method.

Looking at what λ does in terms of singular values makes it clear. Each singular value normally gets inverted as 1/σ, but with damping added, it becomes σ/(σ² + λ²). When σ is much larger than λ, this is nearly the same as 1/σ; when σ goes to zero, this also goes to zero. It selectively suppresses the amplification along small-singular-value directions without disturbing the rest.

import numpy as np

L1, L2 = 0.20, 0.15


def fk(q):
    return np.array([L1 * np.cos(q[0]) + L2 * np.cos(q[0] + q[1]),
                     L1 * np.sin(q[0]) + L2 * np.sin(q[0] + q[1])])


def jacobian(q):
    s1, c1 = np.sin(q[0]), np.cos(q[0])
    s12, c12 = np.sin(q[0] + q[1]), np.cos(q[0] + q[1])
    return np.array([[-L1 * s1 - L2 * s12, -L2 * s12],
                     [ L1 * c1 + L2 * c12,  L2 * c12]])


J = jacobian(np.radians([30, 1]))          # nearly fully extended, right at the edge of a singularity
dx = np.array([0.001, 0.0])                # want to move 1mm in the x direction

dq = np.linalg.pinv(J) @ dx
print(f"Pseudoinverse  |Δq|={np.linalg.norm(dq):.6f} rad = {np.degrees(np.linalg.norm(dq)):8.4f}°"
      f"  at 100Hz {np.degrees(np.linalg.norm(dq))/0.01:8.1f}°/s"
      f"  actual motion {np.linalg.norm(J @ dq)*1000:.4f}mm")

for lam in (0.001, 0.005, 0.01, 0.05):
    dq = J.T @ np.linalg.solve(J @ J.T + lam * lam * np.eye(2), dx)
    print(f"DLS λ={lam:<6} |Δq|={np.linalg.norm(dq):.6f} rad = {np.degrees(np.linalg.norm(dq)):8.4f}°"
          f"  at 100Hz {np.degrees(np.linalg.norm(dq))/0.01:8.1f}°/s"
          f"  actual motion {np.linalg.norm(J @ dq)*1000:.4f}mm")

Here's the output.

Pseudoinverse  |Δq|=0.626520 rad =  35.8970°  at 100Hz   3589.7°/s  actual motion 1.0000mm
DLS λ=0.001  |Δq|=0.409784 rad =  23.4789°  at 100Hz   2347.9°/s  actual motion 0.7585mm
DLS λ=0.005  |Δq|=0.044071 rad =   2.5251°  at 100Hz    252.5°/s  actual motion 0.5113mm
DLS λ=0.01   |Δq|=0.011702 rad =   0.6705°  at 100Hz     67.0°/s  actual motion 0.5077mm
DLS λ=0.05   |Δq|=0.001394 rad =   0.0799°  at 100Hz      8.0°/s  actual motion 0.4992mm

This table shows exactly what's being traded off.

The pseudoinverse produces the requested 1 millimeter exactly. In exchange, the joints have to move 35.9 degrees.

At λ = 0.01, joint motion drops to 0.67 degrees — a 54x reduction. In exchange, the actual motion is only 0.51 millimeters, about half of what was requested.

Only going halfway looks like a failure, but it isn't. This is the story for one cycle. Next cycle, it looks at the remaining error again and goes about half of that again. The end effector approaches the target a bit slowly, but reliably. The pseudoinverse, by contrast, tries to arrive in a single cycle, pushes the joints past their physical limits, and as a result the end effector ends up somewhere completely wrong instead of anywhere near the target.

There's a practical guideline for choosing λ. If you're using a fixed value, work it out backward from your required precision and maximum joint speed. A better approach is variable damping: λ = 0 far from a singularity, and only ramp λ up once the smallest singular value drops below a threshold. This way you don't lose everyday accuracy while still getting protection near a singularity.

Redundancy and Null Space

One last piece. What happens when there are more joints than task dimensions?

Say we specify only a 2D position with a 3-link planar arm. The Jacobian is 2 rows by 3 columns, with rank 2. By the rank-nullity theorem, the null space's dimension is 3 - 2 = 1.

A joint velocity vector in the null space satisfies J·Δq = 0. The joints move, but the end effector stays put. That's exactly the motion of keeping your hand fixed on a desk while moving your elbow up and down.

The way you use this is a null-space projection.

Δq = J⁺·e + (I - J⁺·J)·z

The first term carries the end effector to the target, and the second term nudges the pose in whatever direction z wants, without affecting the end effector at all. (I - J⁺·J) is the matrix that projects onto the null space.

What you put into z decides how you use the redundancy. Feed in a value that grows the farther a joint sits from the center of its range, and the arm steers away from its limits. Feed in the distance to an obstacle, and the arm dodges its own body while holding the end-effector path fixed. Feed in the manipulability index, and the arm steers itself away from singularities.

import numpy as np

LINKS = (0.20, 0.15, 0.10)


def fk3(q):
    a = np.cumsum(q)
    return np.array([sum(L * np.cos(t) for L, t in zip(LINKS, a)),
                     sum(L * np.sin(t) for L, t in zip(LINKS, a))])


def jacobian3(q):
    a = np.cumsum(q)
    J = np.zeros((2, 3))
    for j in range(3):
        J[0, j] = -sum(LINKS[i] * np.sin(a[i]) for i in range(j, 3))
        J[1, j] =  sum(LINKS[i] * np.cos(a[i]) for i in range(j, 3))
    return J


def solve_ik(target, q0, lam=0.05, max_iter=200, tol=1e-6):
    """Damped least squares iteration. Doesn't diverge even passing through a singularity."""
    q = np.array(q0, dtype=float)
    for k in range(max_iter):
        error = target - fk3(q)
        if np.linalg.norm(error) < tol:
            return q, k
        J = jacobian3(q)
        q = q + J.T @ np.linalg.solve(J @ J.T + lam * lam * np.eye(2), error)
    return q, max_iter


target = np.array([0.30, 0.20])
q, iters = solve_ik(target, [0.1, 0.5, 0.0])
print(f"Converged in {iters} iterations  θ={np.round(np.degrees(q), 4)}°  FK check={np.round(fk3(q), 8)}")

J = jacobian3(q)
print(f"J shape {J.shape}, rank {np.linalg.matrix_rank(J)}, null-space dimension {3 - np.linalg.matrix_rank(J)}")

N = np.eye(3) - np.linalg.pinv(J) @ J     # null-space projection matrix
z = np.array([1.0, 1.0, 1.0])             # try feeding in an arbitrary direction
dq_null = N @ z
print(f"Null-space direction Δq = {np.round(dq_null, 6)}")
print(f"End-effector velocity in this direction J·Δq = {np.round(J @ dq_null, 12)}")

step = 0.05 * dq_null / np.linalg.norm(dq_null)
print(f"After moving the joints by {np.round(np.degrees(step), 4)}°")
print(f"  end-effector displacement = {np.linalg.norm(fk3(q + step) - fk3(q)) * 1000:.6f} mm")

Here's the output.

Converged in 6 iterations  θ=[-7.5005 67.2642 15.0619]°  FK check=[0.29999971 0.19999971]
J shape (2, 3), rank 2, null-space dimension 1
Null-space direction Δq = [ 0.038069 -0.231649  0.463807]
End-effector velocity in this direction J·Δq = [0. 0.]
After moving the joints by [ 0.2098 -1.2766  2.556 ]°
  end-effector displacement = 0.059555 mm

The three joints moved 0.21°, -1.28°, and 2.56° respectively, and the end effector barely moved at all — 0.06 millimeters. The reason it isn't exactly zero is that the null space is the tangent direction at that one point, and a finite step leaves a second-order error. Halve the step and this error drops to a quarter.

It's also worth noting the iteration converged in 6 steps. Despite starting fairly far from the target, it lands within micrometers in just six iterations. Newton-family methods converge this fast. Of course, if it passes through a singularity, or the target lies outside the workspace, the story changes, so real code always needs an iteration-count cap and a final error check. If you noticed that the function above silently returns a value even when it hits max_iter, you noticed correctly — production code needs to let the caller distinguish that case.

Conclusion — Inverse Kinematics Isn't a Problem of Finding an Answer, It's a Problem of Choosing One

One fact runs through this whole post: inverse kinematics has several answers, or none at all.

So most of the actual work isn't solving an equation — it's specifying which answer you want. Elbow-up or elbow-down? Discard solutions beyond joint limits? Pick the solution closest to the previous pose? If there's redundancy, what do you spend the extra freedom on? Near a singularity, do you give up accuracy or joint speed?

Fail to make these choices explicit in code, and the choice doesn't disappear — floating-point arithmetic order picks for you instead. And that choice changing cycle to cycle is exactly why an arm suddenly flips over.

One last thing to remember when using a numerical method: the iteration converging and the answer being usable are two different things. Whether the converged angles fall within joint limits, whether the path along the way passes through the arm's own body, whether that pose's condition number is manageable — all of this needs to be checked separately.

That covers "where should it go." Actually moving smoothly to that angle continues in the control loop post. And how far you need to study the tools that showed up in this post — the pseudoinverse, singular value decomposition, condition number — is laid out in order in the math you need for robotics post.