deepxde icon indicating copy to clipboard operation
deepxde copied to clipboard

Biharmonic equation

Open liangchun0206 opened this issue 2 years ago • 3 comments

Hi @lululxvi. Is it possible to implement the biharmonic equation (u_xxxx+2u_xxyy+u_yyyy=f) numerically by DeepXDE? Here is my code, but it does not work.

from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde from deepxde.backend import tf from scipy.special import gamma

def pde(x, y): dy_xx = dde.grad.hessian(y, x, i=0, j=0) dy_yy = dde.grad.hessian(y, x, i=1, j=1) dy_xxxx= dde.grad.hessian(y, dy_xx, i=0, j=0) dy_xxyy= dde.grad.hessian(y, dy_xx, i=1, j=1) dy_yyyy= dde.grad.hessian(y, dy_yy, i=1, j=1)

f=np.exp(x[:, 0:1])*np.sin(np.pi*x[:, 1:2])-2*np.pi**2*np.sin(np.pi*x[:, 1:2])+np.pi**4*np.exp(x[:, 0:1])*np.sin(np.pi*x[:, 1:2])
return dy_xxxx+2*dy_xxyy+dy_yyyy-f

def boundary_l(x, on_boundary): return on_boundary

def boundary_r1(x, on_boundary): return on_boundary and np.isclose(x[0], 0)

def boundary_r2(x, on_boundary): return on_boundary and np.isclose(x[0], 1)

def boundary_r3(x, on_boundary): return on_boundary and np.isclose(x[1], 0)

def boundary_r4(x, on_boundary): return on_boundary and np.isclose(x[1], 1)

def func(x): return np.exp(x[:, 0:1])np.sin(np.pix[:, 1:2])

geom = dde.geometry.Rectangle([0,0], [1,1])

bc_l = dde.icbc.DirichletBC(geom, func, boundary_l) bc_r1 = dde.icbc.NeumannBC(geom, lambda X: -np.sin(np.piX[:, 1:2]), boundary_r1) bc_r2 = dde.icbc.NeumannBC(geom, lambda X: np.expnp.sin(np.piX[:, 1:2]), boundary_r2) bc_r3 = dde.icbc.NeumannBC(geom, lambda X: -np.pinp.exp(X[:, 0:1]), boundary_r3) bc_r4 = dde.icbc.NeumannBC(geom, lambda X: -np.pi*np.exp(X[:, 0:1]), boundary_r4)

data = dde.data.PDE( geom, pde, [bc_l,bc_r1,bc_r2,bc_r3,bc_r4], num_domain=1200, num_boundary=120, solution=func,num_test=1500) net = dde.maps.FNN([2] + [50] * 4 + [1], "tanh", "Glorot uniform") model = dde.Model(data, net)

model.compile("adam", lr=0.001) model.train(epochs=50000) #model.compile("L-BFGS-B") losshistory, train_state = model.train() dde.saveplot(losshistory, train_state, issave=True, isplot=True)

X = geom.random_points(1000) y_true = func(X) y_pred = model.predict(X) print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred)) np.savetxt("test.dat", np.hstack((X, y_true, y_pred)))

liangchun0206 avatar May 18 '22 17:05 liangchun0206

I think it is doable. But I have found some bug in your code.

  1. The definition of 4th-order derivatives should be something like dy_xxxx=dde.grad.hessian(dy_xx, x, i=0, j=0).
  2. You should use tf instead of np in pde
  3. The func and bc miss some operators. For example, what is np.exp(x[:, 0:1])np.sin(np.pix[:, 1:2])?

forxltk avatar May 18 '22 18:05 forxltk

I have finish the code. Here is it from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde from deepxde.backend import tf from scipy.special import gamma

def pde(x, y): dy_xx = dde.grad.hessian(y, x, i=0, j=0) dy_yy = dde.grad.hessian(y, x, i=1, j=1) dy_xxxx= dde.grad.hessian(dy_xx, x, i=0, j=0) dy_xxyy= dde.grad.hessian(dy_xx, x, i=1, j=1) dy_yyyy= dde.grad.hessian(dy_yy, x, i=1, j=1)

f=tf.exp(x[:, 0:1])*tf.sin(np.pi*x[:, 1:2])-2*np.pi**2*tf.sin(np.pi*x[:, 1:2])+np.pi**4*tf.exp(x[:, 0:1])*tf.sin(np.pi*x[:, 1:2])
return dy_xxxx+2*dy_xxyy+dy_yyyy-f

def boundary_l(x, on_boundary): return on_boundary

def boundary_r1(x, on_boundary): return on_boundary and np.isclose(x[0], 0)

def boundary_r2(x, on_boundary): return on_boundary and np.isclose(x[0], 1)

def boundary_r3(x, on_boundary): return on_boundary and np.isclose(x[1], 0)

def boundary_r4(x, on_boundary): return on_boundary and np.isclose(x[1], 1)

def func(x): return np.exp(x[:, 0:1])np.sin(np.pix[:, 1:2])

geom = dde.geometry.Rectangle([0,0], [1,1])

bc_l = dde.icbc.DirichletBC(geom, func, boundary_l) bc_r1 = dde.icbc.NeumannBC(geom, lambda X: -tf.sin(np.piX[:, 1:2]), boundary_r1) bc_r2 = dde.icbc.NeumannBC(geom, lambda X: np.exp([1])tf.sin(np.piX[:, 1:2]), boundary_r2) bc_r3 = dde.icbc.NeumannBC(geom, lambda X: -np.pitf.exp(X[:, 0:1]), boundary_r3) bc_r4 = dde.icbc.NeumannBC(geom, lambda X: -np.pi*tf.exp(X[:, 0:1]), boundary_r4)

data = dde.data.PDE( geom, pde, [bc_l,bc_r1,bc_r2,bc_r3,bc_r4], num_domain=1200, num_boundary=200, solution=func, num_test=50000) net = dde.maps.FNN([2] + [100] * 2 + [1], "tanh", "Glorot uniform")

model = dde.Model(data, net) model.compile("adam", lr=0.001, metrics=["l2 relative error"]) #model.compile("L-BFGS-B") losshistory, train_state = model.train(epochs=50000) dde.saveplot(losshistory, train_state, issave=True, isplot=True)

X = geom.random_points(5000) y_true = func(X) y_pred = model.predict(X) print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred)) np.savetxt("test.dat", np.hstack((X, y_true, y_pred)))

liangchun0206 avatar May 21 '22 15:05 liangchun0206

is this working or are you getting any errors, is so can you share the error.

NischalKMapakshi avatar Nov 02 '23 07:11 NischalKMapakshi