I’m moving my first steps with NeuralPDE.jl and I came up with this issue. I want to solve a fairly simple PDE:
h(ω) = exp.(-(γ-1)*ω) μω = λ*(ω-ωb) σω = σ eq = 1/100*Dτ(f(τ,ω)) ~ 0.5*σω^2*Dωω(f(τ,ω)) - μω*Dω(f(τ,ω)) - k*f(τ,ω) + h(ω)
where \omega\in[-2, 2] and \tau\in[0,1]. I managed to get a decent solution (just by visual comparison with the finite difference one, see figure below for a cross section at time \tau=1) by using GridTraining:
On the other hand, if I use QuadratureTraining, with say CubatureJLh (one of the suggested algos) I get a small loss function but the solution is nowhere near the real one (same cross section as before):
I tried different cubature methods and played around with some parameters with no luck.
Am I doing something obviously wrong? Is the problem particularly not suited to cubature? Any insight would be appreciated
Here’s some of the code I use to solve the problem if it might help:
τ_min= 0. τ_max = 1.0 τ_span = (τ_min,τ_max) ω_min = -2. ω_max = 2. dω = 0.05 dt_NN = 0.01 dt_BSDE = 0.01 μY = 0.03 σ = 0.03 λ = 0.1 ωb = μY/λ γ = 1.5 k = 0.1 h(ω) = exp.(-(γ-1)*ω) @parameters τ ω @variables f(..) Dω = Differential(ω) Dωω = Differential(ω)^2 Dτ = Differential(τ) μω = λ*(ω-ωb) σω = σ # PDE eq = 1/100*Dτ(f(τ,ω)) ~ 0.5*σω^2*Dωω(f(τ,ω)) - μω*Dω(f(τ,ω)) - k*f(τ,ω) + h(ω) # Initial and boundary conditions bcs = [f(τ_min,ω) ~ h(ω), Dω(f(τ,ω_max)) ~ 0, Dω(f(τ,ω_min)) ~ 0] # Space and time domains domains = [τ ∈ Interval(τ_min,τ_max), ω ∈ Interval(ω_min,ω_max)] τs,ωs = [infimum(d.domain):dω:supremum(d.domain) for d in domains] @named pde_system = PDESystem(eq,bcs,domains,[τ,ω],[f(τ, ω)]) # NN parameters dim =2 hls = dim+50 # Neural network chain = Chain(Dense(dim,hls,Flux.σ), Dense(hls,hls,Flux.σ), Dense(hls,hls,Flux.σ), Dense(hls,1)) strategy = GridTraining([dω,dt_NN]) # strategy = QuadratureTraining(; quadrature_alg = CubatureJLp(), # reltol = 1e-6, abstol = 1e-3, # maxiters = 1_000) discretization = PhysicsInformedNN(chain, strategy) prob = discretize(pde_system,discretization) callback = function (p,l) println("Current loss is: $l") return false end res = Optimization.solve(prob,Flux.ADAM(0.01);callback = callback,maxiters=500) prob = remake(prob,u0=res.u) res = Optimization.solve(prob,Flux.ADAM(0.001);callback = callback,maxiters=2000) prob = remake(prob,u0=res.u) bfgs = OptimizationOptimJL.BFGS() res = Optimization.solve(prob,bfgs ;callback = callback,maxiters=100) prob = remake(prob,u0=res.u) phi = discretization.phi f_NN = [first(Array(phi([τ, ω], res.u))) for τ in τs for ω in ωs] f_NN = reshape(f_NN,(length(ωs),length(τs)))