```python
# https://colab.research.google.com/github/kassbohm/wb-snippets/blob/master/ipynb/HTM_03/Uebungen/3.1.ipynb

# Rounding:
import decimal
from decimal import Decimal as DX
from copy import deepcopy
def iso_round(obj, pv, rounding=decimal.ROUND_HALF_EVEN):
    import sympy
    """
    Rounding acc. to DIN EN ISO 80000-1:2013-08
    place value = Rundestellenwert
    """
    assert pv in set([
        # place value   #  round to:
        1,              #  1
        0.1,            #  1st digit after decimal
        0.01,           #  2nd
        0.001,          #  3rd
        0.0001,         #  4th
        0.00001,        #  5th
        0.000001,       #  6th
        0.0000001,      #  7th
        0.00000001,     #  8th
        0.000000001,    #  9th
        0.0000000001,   # 10th
        ])
    objc = deepcopy(obj)
    try:
        tmp = DX(str(float(objc)))
        objc = tmp.quantize(DX(str(pv)), rounding=rounding)
    except:
        for i in range(len(objc)):
            tmp = DX(str(float(objc[i])))
            objc[i] = tmp.quantize(DX(str(pv)), rounding=rounding)
    return objc

from sympy import *
from sympy.physics.units import *
from mpmath import radians, degrees, pi

l1, l2, l3, b = 1*m, 2*m, 2*m, 2*m
p1 = radians(30)

p2, p3 = var("phi2, phi3")
dp2, dp3 = var("delta_2, delta_3")

c1, s1 = cos(p1), sin(p1)
c2, s2 = cos(p2), sin(p2)
c3, s3 = cos(p3), sin(p3)

xdev = ( l1*c1 + l2*c2 - l3*c3 - b ) / m
ydev = ( l1*s1 + l2*s2 - l3*s3 ) / m

J11, J12 = diff(xdev, p2), diff(xdev, p3)
J21, J22 = diff(ydev, p2), diff(ydev, p3)

# Jacobian:
J = Matrix([[J11, J12],[J21, J22]])

pprint("\n--- Numerical Solution using SymPy's nsolve ---")

# Initial Values:
p2i, p3i = radians(45), radians(90)

sol = nsolve( [xdev, ydev], [p2, p3], [p2i , p3i], dict=True )
sol = sol[0]
sol =  Matrix([sol[p2], sol[p3]])
pprint("\n(phi2, phi3) / deg:")
tmp = degrees(sol)
pprint(iso_round(tmp,0.01))

pprint("\n--- Step-by-Step Newton-Solution ---")
pprint("\n1: Jacobian:")
pprint(J)

# n = 1 Iteration:
pprint("\n2: Init:")
p = Matrix([p2i, p3i])
pprint(iso_round(p, 0.01))

F = Matrix([xdev, ydev])

sub_list = [
    (p2, p2i),
    (p3, p3i),
    ]

pprint("\n3: Jacobian for (phi2, phi3):")
Jn = J.subs(sub_list)
# pprint(Jn)
pprint(iso_round(Jn, 0.01))

pprint("\n4: Deviation Phi for (phi2, phi3):")
Fn = F.subs(sub_list)
pprint(iso_round(Fn, 0.01))

# Newton-Step:
pprint("\n5: New (phi2, phi3):")
dp = Matrix([dp2, dp3])
eq = Eq(Jn*dp, -Fn)
sol = solve(eq,[dp2, dp3], dict=True)
sol = sol[0]
dp = Matrix([sol[dp2], sol[dp3]])
p += dp

pprint("\nIn rad (radians):")
pprint(iso_round(p, 0.01))

pprint("\nIn deg (degrees):")
pprint(iso_round(degrees(p), 0.01))

# --- Numerical Solution using SymPy's nsolve ---
#
# (phi2, phi3) / deg:
# ⎡48.16⎤
# ⎢     ⎥
# ⎣84.25⎦
#
# --- Step-by-Step Newton-Solution ---
#
# 1: Jacobian:
# ⎡-2.0⋅sin(φ₂)  2.0⋅sin(φ₃) ⎤
# ⎢                          ⎥
# ⎣2.0⋅cos(φ₂)   -2.0⋅cos(φ₃)⎦
#
# 2: Init:
# ⎡0.79⎤
# ⎢    ⎥
# ⎣1.57⎦
#
# 3: Jacobian for (phi2, phi3):
# ⎡-1.41  2.0⎤
# ⎢          ⎥
# ⎣1.41   0.0⎦
#
# 4: Deviation Phi for (phi2, phi3):
# ⎡0.28 ⎤
# ⎢     ⎥
# ⎣-0.09⎦
#
# 5: New (phi2, phi3):
#
# In rad (radians):
# ⎡0.85⎤
# ⎢    ⎥
# ⎣1.47⎦
#
# In deg (degrees):
# ⎡48.48⎤
# ⎢     ⎥
# ⎣84.43⎦

```

9c5e1ce3892c9886087f72e36a50e84e4dca3634
ipynb
Jupyter Notebook
ipynb/HTM_03/Uebungen/3.1.ipynb
kassbohm/wb-snippets
f1ac5194e9f60a9260d096ba5ed1ce40b844a3fe
MIT
ipynb/HTM_03/Uebungen/3.1.ipynb
kassbohm/wb-snippets
f1ac5194e9f60a9260d096ba5ed1ce40b844a3fe
MIT
ipynb/HTM_03/Uebungen/3.1.ipynb
kassbohm/wb-snippets
f1ac5194e9f60a9260d096ba5ed1ce40b844a3fe
MIT
Qwen/Qwen-72B
1. YES
2. YES