神经网络近似解析结果

神经网络近似解析结果.png
(2+1)模板来这里

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux, Plots, Quadrature, Cubature,DiffEqFlux        
import ModelingToolkit: Interval, infimum, supremum
CUDA.device()

输入基本信息

@parameters t, x
@variables u(..)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Dt = Differential(t)

#2D PDE
C=1
eq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))

# Initial and boundary conditions
bcs = [u(t,0) ~ 0.,# for all t > 0
       u(t,1) ~ 0.,# for all t > 0
       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]
# Discretization
dx = 0.1
# Neural network
chain = FastChain(FastDense(2,3,tanh),FastDense(3,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
#Chain(Dense(1, 64, σ), Dense(64, 64, σ) , Dense(5, 2))

生成数据点

discretization = PhysicsInformedNN(chain, GridTraining(dx); init_params = initθ)

生成最优化问题 即 LOSS问题

@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
prob = discretize(pde_system,discretization)

求解最优化问题

#每次输出loss值
cb = function (p,l)
    println("Current loss is: $l")
    return false
end

# optimizer
res = GalacticOptim.solve(prob,Optim.BFGS(); cb = cb, maxiters=10)
#,use_gpu = true

• 搞出来解析解

phi = discretization.phi
first(phi([t,x],res.minimizer)) #(phi([t,x],res.minimizer))[1]

绘图

ts,xs = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
analytic_sol_func(t,x) =  sum([(8/(k^3*pi^3)) * sin(k*pi*x)*cos(C*k*pi*t) for k in 1:2:50000])

u_predict = reshape([first(phi([t,x],res.minimizer)) for t in ts for x in xs],(length(ts),length(xs)))
u_real = reshape([analytic_sol_func(t,x) for t in ts for x in xs], (length(ts),length(xs)))

diff_u = abs.(u_predict .- u_real)
p1 = plot(ts, xs, u_real, linetype=:contourf,title = "analytic");
p2 =plot(ts, xs, u_predict, linetype=:contourf,title = "predict");
p3 = plot(ts, xs, diff_u,linetype=:contourf,title = "error");

plot(p1,p2,p3)
p1= plot(xs, u_predict[3],title = "t = 0.1");
p2= plot(xs, u_predict[11],title = "t = 0.5");
p3= plot(xs, u_predict[end],title = "t = 1");
plot(p1,p2,p3)