Matrix form of "A simple ODE system"

[louisricoeur]

Dear @lululxvi and team,

I have been trying to replicate the "A simple ODE system" example with the ODE specified in a matrix form - before adapting it to my specific problem. I got the message: InvalidArgumentError: Incompatible shapes: [39,1] vs. [2,39] [[{{node sub_200}}]]

The only part I changed was the one involving the ODE system:

A = np.array([[0, 1], [-1, 0]])
A_tf = tf.convert_to_tensor(A, dtype=tf.float32)

def ode_deepxde(x, y):
    lhs = dde.grad.jacobian(y, x)
    rhs = tf.matmul(A_tf, y, transpose_b=True)
    return lhs - rhs

Now A_tf is a 2x2 matrix, y should be a 1x2 vector and dde.grad.jacobian(y, x) should be a 2x1 vector. Could you please help me trouble shoot this issue? Thank you.

I have already looked the Lotka-Volterra example where tf.matmul is used, but could not adapt it successfully.

[haison19952013]

Instead of using the matrix form, I think you just need to derive that 2x2 matrix into 2 equations.It means you will have 2 residuals. Plz following the tutorial.

[louisricoeur]

Thank you for your answer.

However, the objective is then to apply the technique to my specific application, where I have a linear system with 8x8 matrix which is the outcome of an algorithm. So it would be way more convenient to be able to use the matrix directly, rather than manually specify dy1_x etc.

So your suggestion does not solve the problem unfortunately.

[praksharma]

Please put the full trackback of the error inside the code snippet i.e. <> icon in the toolbar.

When you say "ODE specified in a matrix form" i guess you are talking about the global equilibrium equation for example, stiffness*displacement-force = 0? Is that correct?

[louisricoeur]

Dear @praksharma, thank you for coming back to me regarding this topic.

• When I say "ODE specified in a matrix form", I mean any ODE written in the form dy/dt = Ay, where y is a nx1 column vector and A is a mxn matrix.

For the specific system given in "A simple ODE system", it means writing the ODE as dy/dt = Ay where: A = \begin{matrix} 0 & 1\ -1 & 0 \end{matrix}.

• Regarding the full trackback of the error, please find it below:

InvalidArgumentError Traceback (most recent call last) File /usr/local/anaconda3/envs/DeepXDE/lib/python3.10/site-packages/tensorflow/python/client/session.py:1377, in BaseSession._do_call(self, fn, *args) 1376 try: -> 1377 return fn(*args) 1378 except errors.OpError as e:

File /usr/local/anaconda3/envs/DeepXDE/lib/python3.10/site-packages/tensorflow/python/client/session.py:1360, in BaseSession._do_run.<locals>._run_fn(feed_dict, fetch_list, target_list, options, run_metadata) 1359 self._extend_graph() -> 1360 return self._call_tf_sessionrun(options, feed_dict, fetch_list, 1361 target_list, run_metadata)

File /usr/local/anaconda3/envs/DeepXDE/lib/python3.10/site-packages/tensorflow/python/client/session.py:1453, in BaseSession._call_tf_sessionrun(self, options, feed_dict, fetch_list, target_list, run_metadata) 1451 def _call_tf_sessionrun(self, options, feed_dict, fetch_list, target_list, 1452 run_metadata): -> 1453 return tf_session.TF_SessionRun_wrapper(self._session, options, feed_dict, 1454 fetch_list, target_list, 1455 run_metadata)

InvalidArgumentError: Incompatible shapes: [39,1] vs. [2,39] [[{{node sub}}]]

` During handling of the above exception, another exception occurred:

InvalidArgumentError Traceback (most recent call last) ... File "/usr/local/anaconda3/envs/DeepXDE/lib/python3.10/site-packages/tensorflow/python/framework/ops.py", line 3754, in _create_op_internal ret = Operation( File "/usr/local/anaconda3/envs/DeepXDE/lib/python3.10/site-packages/tensorflow/python/framework/ops.py", line 2133, in init self._traceback = tf_stack.extract_stack_for_node(self._c_op) `

• For the record, the full code I am using is:

import numpy as np
A = np.array([[0, 1], [-1, 0]])

def ode_deepxde(x, y):
    lhs = dde.grad.jacobian(y, x)
    rhs = tf.matmul(A_tf, y, transpose_b=True)
    return lhs - rhs

geom = dde.geometry.TimeDomain(0, Tf)
def boundary(_, on_initial):
    return on_initial
ic1 = dde.icbc.IC(geom, lambda x: 0, boundary, component=0)
ic2 = dde.icbc.IC(geom, lambda x: 1, boundary, component=1)
data = dde.data.PDE(geom, ode_deepxde, [ic1, ic2], 35, 2, num_test=100)
layer_size = [1] + [50] * 3 + [A.shape[0]]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=0.001)
losshistory, train_state = model.train(iterations=20000, display_every=100)

[praksharma]

InvalidArgumentError: Incompatible shapes: [39,1] vs. [2,39].

You might want to replace

rhs = tf.matmul(A_tf, y, transpose_b=False)

I didn't run your code.

[louisricoeur]

Just implementing the change you suggested leads to another, similar error:

InvalidArgumentError: Matrix size-incompatible: In[0]: [2,2], In[1]: [39,2] [[{{node MatMul}}]]

Changing this same line in the following way does not lead to an error message, but the solution is completely off, which should not be the case given that writting the ODE in a matrix form is (mathematically) equivalent to specifying it line by line and that the same code works when the ODE is specified line by line:

rhs = tf.matmul(y, A_tf, transpose_b=True)

Same comment goes when changing this same line to:

rhs = tf.matmul(y, A_tf, transpose_b=False)

I am a bit puzzled at what to try next.

[lululxvi]

Note that y has the shape batch_size x N.

[riccardotomada]

I just went through the same problem. If it may help anyone, here it is the full code which can be easily extended to the case in which N >> 1.

import deepxde as dde
import numpy as np

A_tf = np.array([[0, 1], [-1, 0]], dtype = 'float64')
n = 2

def ode_deepxde(x, y):
    lhs0 = dde.grad.jacobian(y, x, i = 0, j = 0)
    for a in range(1,n):
        lhsa = dde.grad.jacobian(y, x, i = a, j = 0)
        lhs  = tf.squeeze(tf.stack((lhs0, lhsa), axis = 1), axis = 2)
    rhs = tf.einsum('jk,kl->jl', A_tf, tf.transpose(y))
    return lhs - tf.transpose(rhs)

geom = dde.geometry.TimeDomain(0, 10)

def boundary(x, on_initial):
    return np.isclose(x[0], 0)

ic1 = dde.icbc.IC(geom, lambda x: 0, boundary, component=0)
ic2 = dde.icbc.IC(geom, lambda x: 1, boundary, component=1)

def func(x):
    return np.hstack((np.sin(x), np.cos(x)))

data = dde.data.PDE(geom, ode_deepxde, [ic1, ic2], 35, 2, num_test=100, solution=func)

layer_size = [1] + [50] * 3 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr=0.001, metrics=["l2 relative error"])
losshistory, train_state = model.train(iterations=20000, display_every=100)