# Parametrization of an ODE coefficient

Hi, I’m trying to solve a system of two coupled equations. The system (the Hopf oscillator), has two coefficients: a and w. I have succesfully trained a model capable of solve the system for a=0.6 and w=1, but now I want to parameterize it with respect to a in the range (-1,1).
Following a tutorial (Parameterized 3D Heat Sink - NVIDIA Docs), I’ve modified my scripts to add the parameterization but it seem like I’m not getting the desired output.
Can someone tell me what I’m doing wrong? And also, how can I plot the results?

``````class Hopf(PDE):
name = "Hopf"

def __init__(self, a=0.6, w=1):
self.a = a
self.a = w

t = Symbol("t")
input_variables = {"t": t}

x = Function("x")(*input_variables)
y = Function("y")(*input_variables)

if type(a) is str:
a = Function(a)(*input_variables)
elif type(a) in [float, int]:
a = Number(a)
if type(w) is str:
w = Function(w)(*input_variables)
elif type(w) in [float, int]:
w = Number(w)

self.equations = {}
self.equations["ode_x"] = x.diff(t)-(a*x - w*y - x*(x**2 + y**2))
self.equations["ode_y"] = y.diff(t)-(a*y + w*x - y*(x**2 + y**2))
``````
``````@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
# make list of nodes to unroll graph on
a, w = 0.6, 1
hopf = Hopf(a=a, w=w)
hopf_net = instantiate_arch(
input_keys=[Key("t"), Key("a")],
output_keys=[Key("x"), Key("y")],
cfg=cfg.arch.fully_connected,
activation_fn=Activation.SIN
)
nodes = hopf.make_nodes() + [hopf_net.make_node(name="hopf_network")]

# make geometry
geo = Point1D(0)
t_min = 0.0
t_max = 20.0
t_symbol = Symbol("t")
a_symbol = Symbol("a")
x = Symbol("x")
y = Symbol("y")
time_range = {t_symbol: (t_min, t_max)}

ic_param = {t_symbol: 0, a_symbol: (-1.0, 1.0)}
initial_param = {t_symbol: (t_min, t_max), a_symbol: (-1.0, 1.0)}

# make domain
domain = Domain()

# initial conditions
IC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"x": 1.0, "y": 0},
batch_size=cfg.batch_size.IC,
# parameterization={t_symbol: 0},
parameterization=ic_param,
)

# solve over given time period
interior = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"ode_x": 0.0, "ode_y": 0.0},
batch_size=cfg.batch_size.interior,
# parameterization=time_range,
parameterization=initial_param,
)

# make solver
slv = Solver(cfg, domain)

# start solver
slv.solve()
``````

Thank you!

Thanks for trying out Modulus. Glad you’ve got the non-parameterized version working! Parameterized solutions can be a lot more tricky but much more useful in the end.

I’ve modified my scripts to add the parameterization but it seem like I’m not getting the desired output. Can someone tell me what I’m doing wrong?

Can you elaborate by what you mean here? For example, Is the convergence bad? Is it unstable? Just not running?

And also, how can I plot the results?

Add a inferencer or validator (for example in the annular ring problem). Additional post processing information in the Modulus-Sym docs!

I wasn’t sure what the problem was, but then I realized that the parameterization wasn’t doing anything (the model was still being trained for one value of ‘a’).

But luckily, I was able to fix the code, so now everything is working as intended. Below I show the full code in case anyone finds it useful.

``````class Hopf(PDE):
name = "Hopf"

def __init__(self, w=1):
self.a = w

t = Symbol("t")
a = Symbol("a")
input_variables = {"t": t, "a": a}

x = Function("x")(*input_variables)
y = Function("y")(*input_variables)

if type(w) is str:
w = Function(w)(*input_variables)
elif type(w) in [float, int]:
w = Number(w)

self.equations = {}
self.equations["ode_x"] = x.diff(t)-(a*x - w*y - x*(x**2 + y**2))
self.equations["ode_y"] = y.diff(t)-(a*y + w*x - y*(x**2 + y**2))
``````
``````@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:

# instantiate ode and network
hopf = Hopf()
hopf_net = instantiate_arch(
input_keys=[Key("t"), Key("a")],
output_keys=[Key("x"), Key("y")],
cfg=cfg.arch.fully_connected,
activation_fn=Activation.SIN
)

# make list of nodes to unroll graph on
nodes = hopf.make_nodes() + [hopf_net.make_node(name="hopf_network")]

# make geometry
geo = Point1D(0)

# create variables
x, y = Symbol("x"), Symbol("y")

# create parameters
a = Symbol("a")
a_range = (cfg.custom.a_min, cfg.custom.a_max)
a_samples = cfg.custom.a_samples

# define time
t_min, t_max = cfg.custom.t_min, cfg.custom.t_max
t = Symbol("t")
t_range = (t_min, t_max)

# initial conditions
x_0 = 1.0
y_0 = 0.0

# make ranges
ic_param_ranges = {t: 0, a: a_range}
interior_param_ranges = {t: t_range, a: a_range}

# make domain
domain = Domain()

# initial conditions
IC = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"x": x_0, "y": y_0},
batch_size=cfg.batch_size.IC,
parameterization=Parameterization(ic_param_ranges),
)

# solve over given time period
interior = PointwiseBoundaryConstraint(
nodes=nodes,
geometry=geo,
outvar={"ode_x": 0.0, "ode_y": 0.0},
batch_size=cfg.batch_size.interior,
parameterization=Parameterization(interior_param_ranges),
)

# for scipy solution
def ode(t, vars, a=0.6, w=1):
x, y = vars
dxdt = a*x - w*y - x*(x**2 + y**2)
dydt = a*y + w*x - y*(x**2 + y**2)
return [dxdt, dydt]
deltaT = 0.001
t_scipy = np.arange(t_min, t_max, deltaT)

# create data for each param value
values_a = np.linspace(*a_range, a_samples)
for sample_a in values_a:
sample_a = float(sample_a)
specific_param_ranges = {
a: sample_a,
}

# add metrics for front pressure
geo_param_ranges = {
**specific_param_ranges,
}

sol = solve_ivp(partial(ode, a=sample_a), (t_min, t_max), [x_0, y_0], t_eval=t_scipy)
invar_numpy = {
"t": np.expand_dims(t_scipy, axis=-1),
"a": np.full_like(np.expand_dims(t_scipy, axis=-1), sample_a)
}
outvar_numpy = {
"x": np.expand_dims(sol.y[0].T, axis=-1),
"y": np.expand_dims(sol.y[1].T, axis=-1),
}

validator = PointwiseValidator(
nodes=nodes, invar=invar_numpy, true_outvar=outvar_numpy, batch_size=1024
)

print(sample_a)

# make solver
slv = Solver(cfg, domain)

# start solver
slv.solve()

if __name__ == "__main__":
run()
``````

I think now I’m having another problem.

I trained a model for ‘a’ in the range [-0.1, 0.1], and got a small loss (≈2e-5). But the results differ a lot from the reference solutions.

For a=0.1:

Is there something that I’m still doing wrong?

This looks like a convergence problem. Its likely the (wrong) prediction does satisfy the PDE loss in some loose form so the error appears small. Consider trying weighting your ODE loss constraint more via the lambda weighting parameter, modifying your network, train more on the earlier time steps than the later or some of the tricks we mention in our user guide.

Personally I would start training a smaller time range (say 0 - 10) and see if it can converge there. If so you know your loss is working, its a matter of how you train it for longer roll outs.

Thank you a lot for the help!! I tried what you suggested: I got the same loss (≈2e-5), but now the results are very good. For a=0.1:

Any tips regarding how to train the model for longer? Because the results are already starting to get worse when I train it for the range 0-20.

Long time predictions is where a lot of experimentation is often required.

• You may have to slowly unroll to longer time scales
• Consider trying some of the other models we support such as the Fourier Networks.
• Are you scaling your time range? Typically a range from 0 to 100 would lead to problems for a neural network.
• May need to train for a longer time with more training points.
• Etc.

Also have a browse of the Recommended Practices section of the docs for some tricks to consider. Good luck!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.