Issue with Solving Second-Order Differential Equation Using DeepXDE

[power817]

Hello,

I am currently working on solving a second-order differential equation using DeepXDE, but I'm encountering some difficulties. The differential equation I'm dealing with is specified as,


and I've set the initial condition as,

I've written some code to tackle this problem, but I'm running into an issue where the second loss always converges to 1, and the solution p(r) I'm getting does not seem to be correct. I'm not sure where the problem lies in my approach or code, and I would greatly appreciate any insights or advice on how to resolve this.

Additionally, I'm unsure how to effectively plot or check the values of the solution (p(r)) after obtaining it. If anyone could guide me on how to visualize the solution or verify its accuracy, that would be incredibly helpful.

Here is the code snippet I'm working with:

def dp(r, p): return dde.grad.jacobian(p, r)

def ddp(r, p): return dde.grad.hessian(p, r)

def At(r): return 0.296977 - 0.296977/r

def U(r): return -2.96977 + 0.0220488/(r**2) + 1.94772/r + r**2

def dU(r): return -0.0440976/(r**3) - 1.94772/(r**2) + 2*r

def M(r): return -1 + 0.0116676/(r**8) - 0.264586/(r**4)

def pde(r, p): dp_r = dp(r, p) dp_rr = ddp(r, p) return (-At(r)**2/(U(r)**2) + M(r)/U(r))*p - (2/r + dU(r)/U(r))*dp_r - dp_rr

def boundary_left(r, on_boundary): return on_boundary and np.isclose(r[0], 1.001)

def boundary_right(r, on_boundary): return on_boundary and np.isclose(r[0], 10000)

def b1(r, p): return p - 1

def b2(r, p): dp_x = dp(r, p) return 122.204 * M(-0.17639) - dp_x

geom = dde.geometry.Interval(1.001, 10000)

bc1 = dde.OperatorBC(geom, lambda r, p, _: b1(r, p), boundary_left)
bc2 = dde.OperatorBC(geom, lambda r, p, _: b2(r, p), boundary_left)

data = dde.data.PDE(geom, pde, [bc1, bc2], num_domain = 1000, num_boundary = 1, num_test = 100)

layer_size = [1] + [20] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"

net = dde.maps.FNN(layer_size, activation, initializer)

model = dde.Model(data, net)
model.compile("adam", lr = 0.001)

losshistory, train_state = model.train(epochs = 10000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
`

I'm looking forward to any suggestions or advice you might have. Thank you in advance for your help!