Commit 317d82bc authored by Frank Hellmann's avatar Frank Hellmann
Browse files

second root finding thing added

parent ea58dd96
......@@ -124,8 +124,24 @@ def define_root_rhs(total_length, rhs):
return np.concatenate((domega, a.view(np.float64)))
return root_rhs
# Another try at a reasonable root finding function. This one just adds omega_global as a variable,
# and adds a zero to the return
def define_root_rhs_omega(total_length, rhs):
zero = np.zeros(1, dtype=np.float64)
def root_rhs(y):
dydt = rhs(y[:-1], 0.)
v = y[:total_length].view(np.complex128)
dv = dydt[:total_length].view(np.complex128)
domega = dydt[total_length:]
omega_global = y[-1]
a = 1.j * omega_global * v - dv
return np.concatenate((zero, domega, a.view(np.float64)))
return root_rhs
def define_network_rhs_codegen(node_list, Y):
network_rhs_numba = None
assert len(node_list) == len(Y)
length = len(Y)
total_length = 2*length
......@@ -189,23 +205,25 @@ if __name__ == "__main__":
root_rhs_1 = lambda x : rhs(x, 0.)
root_rhs = define_root_rhs(2*len(Y), rhs)
root_rhs = define_root_rhs_omega(2*len(Y), rhs)
ic = np.ones(6, dtype=np.float64) * 0.5
guess = np.ones(7, dtype=np.float64) * 0.5
guess[-1] = 0.
print(root_rhs(ic))
print(root_rhs(guess))
from scipy.optimize import root
result = root(root_rhs, ic)
result = root(root_rhs, guess)
if result.success:
print(rhs(result.x, 0.))
ic = result.x
print(root_rhs(ic))
print(rhs(result.x[:-1], 0.))
ic = result.x[:-1]
print(root_rhs(result.x))
# ic[5] += 1.
else:
print("failed")
ic = guess[:-1]
times = np.linspace(0., 100., 1000)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment