Решение нелинейной системы уравнений с комплексными корнями

Question

Я хочу решить нелинейную систему уравнений:

система уравнений

И у меня есть u и v. Вот мой код:

import numpy as np

импортировать matplotlib.pyplot как plt

Dt = 0.002              # timestep Delta t
w = 1
w2 = w
r = 0.33
Tk = 2 * np.pi *r / w
T = Tk

ge = 0
phi = 0.5 / T
q0 = 1
p0 = 1
u0 = np.sqrt(1/2)
v0 = np.sqrt(1/2)

t_start = 0             # starttime
t_end = 40          # endtime
NT = 200
n_steps = int(round((t_end-t_start)/Dt))    # number of timesteps
print(n_steps)
Dt = T / (n_steps + 1)
print(Dt)

q_arr = np.zeros(n_steps + 1, dtype=complex)   # create an array of zeros for q
p_arr = np.zeros(n_steps + 1, dtype=complex)    # create an array of zeros for p
u_arr = np.zeros(n_steps + 1, dtype=complex)   # create an array of zeros for u
v_arr = np.zeros(n_steps + 1, dtype=complex)    # create an array of zeros for v
t_arr = np.zeros(n_steps + 1)   # create an array of zeros for t
t_arr[0] = t_start              # add starttime to array
q_arr[0] = q0              # add initial value of q to array
p_arr[0] = p0              # add initial value of p to array
u_arr[0] = u0              # add initial value of u to array
v_arr[0] = v0              # add initial value of v to array
Q, P, U, V = [q0], [p0], [u0], [v0]

# Euler's method
for j in range(NT):
    
    if j != 0:
        q_arr[0] = q_arr[-1]
        p_arr[0] = p_arr[-1]
        u_arr[0] = u_arr[-1]
        v_arr[0] = v_arr[-1]
        
    for i in range (1, n_steps + 1):
        
        q = q_arr[i-1]
        p = p_arr[i-1]
        u = np.cos(phi) * u_arr[i - 1] - 1j * np.sin(phi) * v_arr[i-1]
        v = -1j * np.sin(phi) * u_arr[i-1] + np.cos(phi) * v_arr[i-1] 
        t = t_arr[i-1]

        dpdt = -w**2 * q - ge * (np.conj(u) * v + u * np.conj(v))  # calculate the derivative of p
        dqdt = p                                                   # calculate the derivative of q
        dudt = -1j * ((w2 / 2) * u + ge * q * v)                   # calculate the derivative of u
        dvdt = -1j * (-(w2 / 2) * v + ge * q * u)                  # calculate the derivative of v
        

        q_arr[i] = q + Dt*dqdt  # calc. q at next timestep,add to array
        p_arr[i] = p + Dt*dpdt  # calc. p at next timestep,add to array
        u_arr[i] = u + Dt*dudt  # calc. u at next timestep,add to array
        v_arr[i] = v + Dt*dvdt  # calc. v at next timestep,add to array
        t_arr[i] = t + Dt       # add new value of t to array
    Q.append(q_arr[-1])
    P.append(p_arr[-1])
    U.append(u_arr[-1])
    V.append(v_arr[-1])
#print(np.real(u_arr))

Я попытался аналитически решить \theta и \phi:

аналитическое решение для тета и фи.

Далее я хочу получить массивы тета и фи:

thett = array([complex(asin(2 * u * v)) for u, v in zip(U, V)])
phii = array([complex(-2j * ln(u / cos(asin(2 * u * v)))) for u, v in zip(U, V)])

Но я обнаружил, что есть ошибка в решении для u и v. В системе Mathematica решение для y и at имеет вид:

математическое решение

и питон дал мне:

решение на питоне.

Я пытаюсь понять, в чем проблема.