Ashutosh Bharambe / Aug 23 2021
GSoC 2021 NumFOCUS : Using the Numeric Integral and solving simple equations
The last blog covers how the Symbolic Integral was added to Julia Symbolics ecosystem, this blog covers the numeric evaluation of the integral and making it compatible with NeuralPDE.jl to solve Integro-Differential Equations.
We start of by adding the necessary packages. Note : I am adding the package NeuralPDE.jl (previously NeuralNetDiffEq.jl) from my forked repository branch, because the pull request is in progress.
] add NeuralPDEusing PkgPkg.add("Flux")Pkg.add("DiffEqFlux")Pkg.add("ModelingToolkit")Pkg.add("DiffEqBase")Pkg.add("GalacticOptim")Pkg.add("Optim")Pkg.add("Quadrature")Pkg.add("SciMLBase")Pkg.add("DomainSets")using Flux, DiffEqFlux, ModelingToolkitusing DiffEqBase, SciMLBaseusing GalacticOptim, Optimusing Quadratureusing DomainSetsusing Plotsusing NeuralPDEWe would start with a simple equation
formula not implementedwith boundary conditions
inline_formula not implemented and inline_formula not implemented
and let's take the domain from x in (0, 2)
We know that the solution of the above equation is u(x) = sin(x)
x u(..)Ix = Integral(x in DomainSets.ClosedInterval(0, pi/2))eq = Ix(u(x)) ~ 1.0bcs = [u(0.0) ~ 0.0, u(1.570) ~ 1.0]domains = [x ∈ Interval(0.0,2.0)]We define the Neural Network chain for the solution of the equation "u"
chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))initθ = Float64.(DiffEqFlux.initial_params(chain))We define our strategy, discretization and the optimizaiton problem
strategy_ = NeuralPDE.GridTraining(0.001)discretization = NeuralPDE.PhysicsInformedNN(chain, strategy_; init_params = nothing, phi = nothing, derivative = nothing, )pde_system = PDESystem(eq,bcs,domains,[x],[u])prob = NeuralPDE.discretize(pde_system,discretization)cb = function (p,l) println("Current loss is: $l") return falseendres = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)Current loss is: 3.2218966211862305
Current loss is: 0.5647337214103231
Current loss is: 0.05065941632700528
Current loss is: 0.01287197469135299
Current loss is: 0.012377104864026081
Current loss is: 0.009810982771725377
Current loss is: 0.008903816662611254
Current loss is: 0.00743997526919167
Current loss is: 0.0007492203397727229
Current loss is: 0.00023809331549511903
Current loss is: 1.7378309452948636e-6
Current loss is: 1.1722953980331601e-8
Current loss is: 6.161841753727519e-12
Current loss is: 2.6348540848573407e-15
Current loss is: 3.0061056792415184e-17
Current loss is: 9.082511224283914e-18
xs = 0.00:0.01:1.5phi = discretization.phiu_predict = [first(phi([x],res.minimizer)) for x in xs]u_real = [sin(x) for x in xs]Flux.mse(u_real, u_predict)plot(xs, u_predict, label = "Predicted Solution")plot!(xs, u_real, label = "Real solution")Let us consider another example with
formula not implementedWe know that the solution would be inline_formula not implemented
Note that we would be using a symbolic upper limit here
Let us write it in the form:
x u(..)Ix = Integral(x in DomainSets.ClosedInterval(0, x))eq = Ix(u(x)) ~ (x^3)/3bcs = [u(0.) ~ 0.0]domains = [x ∈ Interval(0.0,1.00)]chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))initθ = Float64.(DiffEqFlux.initial_params(chain))strategy_ = NeuralPDE.GridTraining(0.05)discretization = NeuralPDE.PhysicsInformedNN(chain, strategy_; init_params = nothing, phi = nothing, derivative = nothing, )pde_system = PDESystem(eq,bcs,domains,[x],[u])prob = NeuralPDE.discretize(pde_system,discretization)res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)xs = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]phi = discretization.phiu_predict = [first(phi([x],res.minimizer)) for x in xs]u_real = [x^2 for x in xs]Current loss is: 0.03693041379451978
Current loss is: 0.013080140445763893
Current loss is: 0.0008532378117857982
Current loss is: 0.0005323377576042182
Current loss is: 0.000532309524927635
Current loss is: 0.0005320024227728011
Current loss is: 0.0005229194137712153
Current loss is: 0.000516389441877494
Current loss is: 0.0005122893463214262
Current loss is: 0.0005109847864859882
Current loss is: 0.0005092952411659875
Current loss is: 0.0005077127812902861
Current loss is: 0.0005072392734359736
Current loss is: 0.0005070771631377893
Current loss is: 0.0005069748305308208
Current loss is: 0.0005069635579767483
Current loss is: 0.0005069606937540017
Current loss is: 0.0005068658796153912
Current loss is: 0.0005037638332184185
Current loss is: 0.0004951562386745019
Current loss is: 0.0004910284938418848
Current loss is: 0.0004482168068829347
Current loss is: 0.00044658861502276675
Current loss is: 0.0004114280628451663
Current loss is: 0.0003523292459125584
Current loss is: 0.0003456835459629048
Current loss is: 0.0003285250768020291
Current loss is: 0.0003042882598483969
Current loss is: 0.0002598032263290685
Current loss is: 0.00017633834901499302
Current loss is: 0.00014282961627793737
Current loss is: 8.95465914774484e-5
Current loss is: 7.770392614452693e-5
Current loss is: 7.281451938150747e-5
Current loss is: 6.182994227391179e-5
Current loss is: 1.0791915633799806e-5
Current loss is: 7.89963227058101e-6
Current loss is: 3.208671269961345e-6
Current loss is: 1.5793561809636117e-6
Current loss is: 7.375476639264053e-7
Current loss is: 6.756176892687944e-7
Current loss is: 6.751597221249109e-7
Current loss is: 6.750881280035624e-7
Current loss is: 6.740865466812316e-7
Current loss is: 4.254988897185772e-7
Current loss is: 3.567767406894336e-7
Current loss is: 3.2666693403430255e-7
Current loss is: 3.1963271785242867e-7
Current loss is: 3.191441936153039e-7
Current loss is: 3.191033255265431e-7
Current loss is: 3.1909859342402837e-7
Current loss is: 3.190872332784509e-7
Current loss is: 3.1908740997022484e-7
Current loss is: 3.1908739997667443e-7
Current loss is: 3.1908740293230564e-7
plot(xs, u_predict, label = "Predicted Solution")plot!(xs, u_real, label = "Real solution")What's Next?
Integro-Differential Equations
Support for multidimensional domains