美文网首页
NeuralPDE:一个求解PDE的包

NeuralPDE:一个求解PDE的包

作者: Neural_PDE | 来源:发表于2021-11-17 12:48 被阅读0次

    Example: Solving 2D Poisson Equation via Physics-Informed Neural Networks

    using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, DiffEqFlux
    using Quadrature, Cubature
    import ModelingToolkit: Interval, infimum, supremum
    
    @parameters x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    
    # 2D PDE
    eq  = Dxx(u(x,y)) + Dyy(u(x,y)) ~ -sin(pi*x)*sin(pi*y)
    
    # Boundary conditions
    bcs = [u(0,y) ~ 0.0, u(1,y) ~ -sin(pi*1)*sin(pi*y),
           u(x,0) ~ 0.0, u(x,1) ~ -sin(pi*x)*sin(pi*1)]
    # Space and time domains
    domains = [x ∈ Interval(0.0,1.0),
               y ∈ Interval(0.0,1.0)]
    # Discretization
    dx = 0.1
    
    # Neural network
    dim = 2 # number of dimensions
    chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
    
    # Initial parameters of Neural network
    initθ = Float64.(DiffEqFlux.initial_params(chain))
    
    discretization = PhysicsInformedNN(chain, QuadratureTraining(),init_params =initθ)
    
    @named pde_system = PDESystem(eq,bcs,domains,[x,y],[u(x, y)])
    prob = discretize(pde_system,discretization)
    
    cb = function (p,l)
        println("Current loss is: $l")
        return false
    end
    
    res = GalacticOptim.solve(prob, ADAM(0.1); cb = cb, maxiters=4000)
    prob = remake(prob,u0=res.minimizer)
    res = GalacticOptim.solve(prob, ADAM(0.01); cb = cb, maxiters=2000)
    phi = discretization.phi
    

    And some analysis:

    xs,ys = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]
    analytic_sol_func(x,y) = (sin(pi*x)*sin(pi*y))/(2pi^2)
    
    u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
    u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
    diff_u = abs.(u_predict .- u_real)
    
    using Plots
    p1 = plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
    p2 = plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
    p3 = plot(xs, ys, diff_u,linetype=:contourf,title = "error");
    plot(p1,p2,p3)
    
    image

    Example: Solving a 100-Dimensional Hamilton-Jacobi-Bellman Equation

    using NeuralPDE
    using Flux
    using DifferentialEquations
    using LinearAlgebra
    d = 100 # number of dimensions
    X0 = fill(0.0f0, d) # initial value of stochastic control process
    tspan = (0.0f0, 1.0f0)
    λ = 1.0f0
    
    g(X) = log(0.5f0 + 0.5f0 * sum(X.^2))
    f(X,u,σᵀ∇u,p,t) = -λ * sum(σᵀ∇u.^2)
    μ_f(X,p,t) = zero(X)  # Vector d x 1 λ
    σ_f(X,p,t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) # Matrix d x d
    prob = TerminalPDEProblem(g, f, μ_f, σ_f, X0, tspan)
    hls = 10 + d # hidden layer size
    opt = Flux.ADAM(0.01)  # optimizer
    # sub-neural network approximating solutions at the desired point
    u0 = Flux.Chain(Dense(d, hls, relu),
                    Dense(hls, hls, relu),
                    Dense(hls, 1))
    # sub-neural network approximating the spatial gradients at time point
    σᵀ∇u = Flux.Chain(Dense(d + 1, hls, relu),
                      Dense(hls, hls, relu),
                      Dense(hls, hls, relu),
                      Dense(hls, d))
    pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
    @time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100,
                                alg=EM(), dt=1.2, pabstol=1f-2)
    

    除此还有很多例子

    相关文章

      网友评论

          本文标题:NeuralPDE:一个求解PDE的包

          本文链接:https://www.haomeiwen.com/subject/kdyltrtx.html