IC for parametric ODE

[npetitclerc]

Hi,
How can I set an initial condition that depends on the input?
I'm trying to add an IC to a parametric ODE, where my parametric y_input = y_output, at t=0.
From the FAQ (#485 (comment)) I thought I could use OperatorBC, but I can't figure out how to use it. Here's what I have tried.
Note: I have 3 parametric inputs and 3 matching outputs.

def boundary_func(input, output, x):
    print(input[:3] )
    print(output)
    return np.abs(input[0:3] - output)

ic1 = dde.icbc.OperatorBC(geom, boundary_func, boundary)

But this fails on the boundary_func on
np.abs(input[0:3] - output)
with:
"ValueError: Dimensions must be equal, but are 4 and 3 for '{{node sub_66}} = Sub[T=DT_FLOAT](strided_slice_121, dense_35/BiasAdd)' with input shapes: [?,4], [?,3]."

I was expecting input[:3] to be a np.array of 3 values, but it's this instead:
Tensor("strided_slice_119:0", shape=(4,), dtype=float32)
and output is
Tensor("strided_slice_120:0", shape=(3,), dtype=float32)

Any idea how to set this IC properly?
And also what are those Tensor variables??

[npetitclerc]

Ok, I figured out 'Tensor' was just a tensorflow Tensor object...
And that converting to numpy array was a bad idea, so using
return input[:, :3] - output
instead doesn't crash, but it gives me significantly worse result, especially near t=0 ... so I'm still doing something wrong.

My first 3 input columns are my 3 parametric IC variables the 4th one is time. So I want to specify that
y_input[:, :3] == y_output, at t=0

My understanding is that 'boundary_func' should be minimized (return values as close to 0 as possible), so I set
return input[:, :3] - output

What I'm I doing wrong here?

[lululxvi]

I am not sure what problem you try to solve.

[npetitclerc]

@lululxvi I'm trying to solve a system of ODEs, to solve the chemical evolution for 3 chemical elements (for this simple test). I have reaction rates at which elements convert into each other.

def ode_system(x, y):
    """ODE system.
    dy1/dx = -1.0 * y1
    dy2/dx = 1.0 * y1 - 0.5* y2
    dy3/dx = 0.5* y2
    """
    y1, y2, y3 = y[:, 0:1], y[:, 1:2], y[:, 2:3]
    dy1_x = dde.grad.jacobian(y, x, i=0)
    dy2_x = dde.grad.jacobian(y, x, i=1)
    dy3_x = dde.grad.jacobian(y, x, i=2)

    return [dy1_x + y1, dy2_x - y1 + 0.5*y2, dy3_x - 0.5*y2]

I can easily solve that with fixed initial conditions, but here I want to pass the initial conditions (initial abundances of my 3 elements) to the network to get the evolved abundances at a different time, aka. make it parametric.
I have been able to do this by setting the geom like this:

geom_t = dde.geometry.TimeDomain(0, 5.0)
geom_p = dde.geometry.Cuboid(xmin=[0.0]*3, xmax=[1.0]*3)
geom = dde.geometry.GeometryXTime(geom_p, geom_t)

and simply ignoring the initial conditions. But to my surprise, I found that the error was largest near t=0. So I wonder if it would help to set my initial conditions as y_input[:, :3] == y_output, at t=0 ? And how to do that.

This is what I tried, but got worse results with:

def boundary_func(input, output, x):
    return tf.abs(input[:,:3] - output)

ic1 = dde.icbc.OperatorBC(geom, boundary_func, boundary)

[haison19952013]

Hi, I think you better use a hard constraint for your outputs. To do that, plz use the apply_output_transform(). For example, suppose you have the inputs $t$, outputs $y_1, y_2$ and different initial conditions $y_1^0, y_2^0$. To make your problem become parametric with $y_1^0, y_2^0$, you have to set your inputs as $[t, y_1^0,y_2^0]$. Then, apply the apply_output_transform() as below:

def modify_output(X, Y):
       t = X[:,0:1]
       y1_0 =  X[:,1:2]
       y2_0 =  X[:,2:3]
       y1 = y[:,0:1]
       y2 = y[:,1:2]
       y1_hard_constraint = (1-t)*y1_0 + t*y1
       y2_hard_constraint = (1-t)*y2_0 + t*y2
       return tf.concat([y1_hard_constraint,y2_hard_constraint ], axis=1)
net.apply_output_transform(modify_output)

Remind that, in my case, $t \in [0, 1]$. However, I think your problem will work with any $t \in [0, \inf]$. My hard constraint scheme here is quite similar to Lagrange interpolation. You can refer in depth about hard constraints in this reference (i.e., section 2.3): https://arxiv.org/pdf/2102.04626.pdf