Ashutosh Bharambe / Aug 23 2021
GSoC 2021 NumFOCUS: Multi Dimensional Integral and Integro Differential Equations
In the last blog we discussed how the single dimension integral is calculated with PINNs, this blog deals with the multi dimensional integral and the integro-differential equations. We'll start by adding the required packages.
using PkgPkg.add("Flux")Pkg.add("DiffEqFlux")Pkg.add("DiffEqBase")Pkg.add("GalacticOptim")Pkg.add("Optim")Pkg.add("Quadrature")Pkg.add("SciMLBase")Pkg.add("DomainSets")Pkg.add("Plots")]add NeuralPDE]dev ModelingToolkitusing Flux, DiffEqFlux, ModelingToolkitusing DiffEqBase, SciMLBaseusing GalacticOptim, Optimusing Quadratureusing DomainSetsWARNING: could not import DistributionsAD._mv_categorical_logpdf into ReverseDiffX
using NeuralPDEMulti Dimensional Integral
Example 1: Definite Integral
Let us start off with a simple definite multi dimensional integral,
inline_formula not implemented
We use our PINN solution inline_formula not implemented to estimate $ 1 - x^2 - y^2 $ Thus, the equation then becomes, inline_formula not implemented
We will provide some boundary conditions for inline_formula not implemented for simplicity:
formula not implementedformula not implementedformula not implementedWe have already discussed how rest of the things are defined, and will see how we define the multi dimensional integral.
x,y u(..)Dx = Differential(x)Dy = Differential(y)Ix = Integral((x,y) in DomainSets.UnitSquare())eq = Ix(u(x,y)) ~ 1/3bcs = [u(0., 0.) ~ 1, Dx(u(x,y)) ~ -2*x , Dy(u(x ,y)) ~ -2*y ]Integral here is defined as Ix = Integral((x,y) in DomainSets.UnitSquare()) . The UnitSquare represents a unit square domain and we can directly use it to define our Integral. The tuple for variables contains variables in the order that they are integrated.
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))initθ = Float64.(DiffEqFlux.initial_params(chain))strategy_ = NeuralPDE.GridTraining(0.1)discretization = NeuralPDE.PhysicsInformedNN(chain, strategy_; init_params = nothing, phi = nothing, derivative = nothing, ) pde_system = PDESystem(eq,bcs,domains,[x,y],[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)xs = 0.00:0.01:1.00ys = 0.00:0.01:1.00phi = discretization.phiu_real = collect(1 - x^2 - y^2 for y in ys, x in xs);u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);┌ Warning: : no method matching get_unit for arguments (SciMLBase.NullParameters,).
└ @ ModelingToolkit /root/.julia/packages/ModelingToolkit/NRQ0J/src/systems/validation.jl:135
Current loss is: 3.5463320215750693
Current loss is: 3.315009641614706
Current loss is: 1.3071218984483521
Current loss is: 1.0832366923639916
Current loss is: 1.0352573604088946
Current loss is: 0.9476226957654448
Current loss is: 0.8960120179288279
Current loss is: 0.882098584779932
Current loss is: 0.8541482141488304
Current loss is: 0.82575633322071
Current loss is: 0.7008364637545456
Current loss is: 0.647617532715176
Current loss is: 0.5387507900730651
Current loss is: 0.45206038709473056
Current loss is: 0.40780736811804325
Current loss is: 0.35184174931018286
Current loss is: 0.322685654032426
Current loss is: 0.25561898826414214
Current loss is: 0.142096747226596
Current loss is: 0.11694088670715294
Current loss is: 0.05087939069224579
Current loss is: 0.04082130960861596
Current loss is: 0.037666745414655575
Current loss is: 0.029215465914965377
Current loss is: 0.02566297096728649
Current loss is: 0.02048262896714069
Current loss is: 0.012435004828455378
Current loss is: 0.008371050643583396
Current loss is: 0.0071056502014713654
Current loss is: 0.00654234932999052
Current loss is: 0.00547525860645941
Current loss is: 0.004506685164656354
Current loss is: 0.004065489456054892
Current loss is: 0.003450364317118752
Current loss is: 0.003208262388937773
Current loss is: 0.003114582310881007
Current loss is: 0.0028105388697045835
Current loss is: 0.002107965347943056
Current loss is: 0.00179947224251424
Current loss is: 0.0016185883912643942
Current loss is: 0.0014439626921880398
Current loss is: 0.0011886183785739138
Current loss is: 0.0010345903946254334
Current loss is: 0.0009407769204129911
Current loss is: 0.0009236767406728457
Current loss is: 0.0008958250894764878
Current loss is: 0.0008802237465781582
Current loss is: 0.0008645784669413438
Current loss is: 0.0008376372064880426
Current loss is: 0.0007920124623028056
Current loss is: 0.0006893556445221664
Current loss is: 0.0005839681997533114
Current loss is: 0.0005344563992558271
Current loss is: 0.0005074382403444835
Current loss is: 0.0004890649377863421
Current loss is: 0.0004757167550136744
Current loss is: 0.0004482318350763225
Current loss is: 0.0004219608874981943
Current loss is: 0.00038107395270317986
Current loss is: 0.00034620858153941274
Current loss is: 0.00025711498736125105
Current loss is: 0.00021857921687027865
Current loss is: 0.00020403786271070166
Current loss is: 0.0001823995273315552
Current loss is: 0.00016798398080083426
Current loss is: 0.00015955583279340528
Current loss is: 0.0001522252800863878
Current loss is: 0.00014191529780499418
Current loss is: 0.00012714736963472217
Current loss is: 0.00010702721881098122
Current loss is: 9.843280382202745e-5
Current loss is: 9.501951838867204e-5
Current loss is: 9.28760268019223e-5
Current loss is: 9.2264590020205e-5
Current loss is: 9.225459071085419e-5
Current loss is: 9.216612836039127e-5
Current loss is: 9.185878339436933e-5
Current loss is: 9.099980191466976e-5
Current loss is: 8.46453794118826e-5
Current loss is: 6.786316847003845e-5
Current loss is: 5.7271003673515025e-5
Current loss is: 5.153878962692276e-5
Current loss is: 4.949847522877105e-5
Current loss is: 4.824618040767939e-5
Current loss is: 4.7441054704634684e-5
Current loss is: 4.693391373873441e-5
Current loss is: 4.6875216427132185e-5
Current loss is: 4.683294633155317e-5
Current loss is: 4.6797841715991556e-5
Current loss is: 4.639147887092654e-5
Current loss is: 4.38063763999099e-5
Current loss is: 3.9065751746976497e-5
Current loss is: 3.478645794119232e-5
Current loss is: 3.357935937481153e-5
Current loss is: 3.1702923492577456e-5
Current loss is: 3.072569747927964e-5
Current loss is: 3.0274619674163255e-5
Current loss is: 3.0077727349701956e-5
Current loss is: 2.9862097330866273e-5
Current loss is: 2.9806019640535074e-5
Current loss is: 2.977332250126614e-5
Now once we are finished with the prediction, lets plot the results
using Plotserror_ = u_predict .- u_realp1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")plot(p1,p2,p3)Example 2: Variable upper Limit
Consider, if we had to use a variable upper limit or lower limit in an integral. Thus, let us look into the integral,
inline_formula not implemented
And we want our PINN solution to estimate inline_formula not implemented
Thus, inline_formula not implemented
We would again add some boundary conditions too.
x,y u(..)Dx = Differential(x)Dy = Differential(y)Ix = Integral((x,y) in DomainSets.ProductDomain(UnitInterval(),ClosedInterval(0 ,x)))eq = Ix(u(x,y)) ~ 5/12bcs = [u(0., 0.) ~ 0, Dy(u(x,y)) ~ 2*y , u(x, 0) ~ x ]Notice here, the Integral is defined as Ix = Integral((x,y) in DomainSets.ProductDomain(UnitInterval(),ClosedInterval(0 ,x))) This is a ProductDomain that we can use to define rectangular and cuboidal domains for our integral.
domains = [x ∈ Interval(0.0,1.00), y ∈ Interval(0.0,1.00)]chain = Chain(Dense(2,15,Flux.σ),Dense(15,1))initθ = Float64.(DiffEqFlux.initial_params(chain))strategy_ = NeuralPDE.GridTraining(0.1)discretization = NeuralPDE.PhysicsInformedNN(chain, strategy_; init_params = nothing, phi = nothing, derivative = nothing, ) pde_system = PDESystem(eq,bcs,domains,[x,y],[u])prob = NeuralPDE.discretize(pde_system,discretization)res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)xs = 0.00:0.01:1.00ys = 0.00:0.01:1.00phi = discretization.phiu_real = collect( x + y^2 for y in ys, x in xs);u_predict = collect(Array(phi([x,y], res.minimizer))[1] for y in ys, x in xs);┌ Warning: : no method matching get_unit for arguments (SciMLBase.NullParameters,).
└ @ ModelingToolkit /root/.julia/packages/ModelingToolkit/NRQ0J/src/systems/validation.jl:135
Current loss is: 4.116698832267988
Current loss is: 1.9300031336141974
Current loss is: 1.0762599571107554
Current loss is: 0.777299749863395
Current loss is: 0.6847827422663673
Current loss is: 0.5371122278980905
Current loss is: 0.5242144781678583
Current loss is: 0.4961384442809944
Current loss is: 0.45325781095415896
Current loss is: 0.4272191431951004
Current loss is: 0.35611895302932073
Current loss is: 0.32858417228315895
Current loss is: 0.2419539231082716
Current loss is: 0.16729640799642634
Current loss is: 0.13078871398417302
Current loss is: 0.1119317149697844
Current loss is: 0.09169758076535031
Current loss is: 0.06908146753242327
Current loss is: 0.06647198114117145
Current loss is: 0.06628618504392889
Current loss is: 0.06600339025970914
Current loss is: 0.06549947568356551
Current loss is: 0.0654478138072065
Current loss is: 0.06539203950779135
Current loss is: 0.06377326964712335
Current loss is: 0.0631273738047732
Current loss is: 0.06254646493484348
Current loss is: 0.06221780042246822
Current loss is: 0.06208052565004679
Current loss is: 0.06199545995774233
Current loss is: 0.06178958631741585
Current loss is: 0.0617137956310438
Current loss is: 0.0616604083465734
Current loss is: 0.06162428361216636
Current loss is: 0.06160002716289341
Current loss is: 0.06158516827083723
Current loss is: 0.06157854381470341
Current loss is: 0.06157392979977729
Current loss is: 0.0615502622048202
Current loss is: 0.061462995126123715
Current loss is: 0.061311204507980895
Current loss is: 0.06102324798911925
Current loss is: 0.06100100606906035
Current loss is: 0.060962687382404186
Current loss is: 0.06094529099877589
Current loss is: 0.06088535529963932
Current loss is: 0.06083019070658782
Current loss is: 0.060824632704265455
Current loss is: 0.060812200967120526
Current loss is: 0.060799964878224276
Current loss is: 0.06079197064443896
Current loss is: 0.060783535830500536
Current loss is: 0.06078007854005499
Current loss is: 0.060773802138916526
Current loss is: 0.06076951065500094
Current loss is: 0.060764233723836186
Current loss is: 0.0607598779366747
Current loss is: 0.06074917676118354
Current loss is: 0.06074144343895895
Current loss is: 0.06071829206718354
Current loss is: 0.06069572980396486
Current loss is: 0.06068979552366958
Current loss is: 0.060677602835786465
Current loss is: 0.060673033573461675
Current loss is: 0.060662750951197886
Current loss is: 0.06065469563181012
Current loss is: 0.06064858798974435
Current loss is: 0.06064515975390796
Current loss is: 0.06064314973189088
Current loss is: 0.060642104313131454
Current loss is: 0.06063886854939005
Current loss is: 0.06063597868393509
Current loss is: 0.060620720926135296
Current loss is: 0.06060990713286811
Current loss is: 0.06059316374927048
Current loss is: 0.060572828541691014
Current loss is: 0.060562899293272514
Current loss is: 0.06055815127040476
Current loss is: 0.06055744763891138
Current loss is: 0.06055718620095305
Current loss is: 0.060557091659755494
Current loss is: 0.06055689859307115
Current loss is: 0.06055656612385959
Current loss is: 0.06055629546317105
Current loss is: 0.060556143727424336
Current loss is: 0.060556097345140385
Current loss is: 0.06055609588188887
Current loss is: 0.06055569513035853
Current loss is: 0.06055498640904425
Current loss is: 0.06055426041207571
Current loss is: 0.060553720029062126
Current loss is: 0.060552059038011
Current loss is: 0.06055117131538573
Current loss is: 0.06054830472814407
Current loss is: 0.06054483203653144
Current loss is: 0.0605356273028702
Current loss is: 0.0605258880879228
Current loss is: 0.06052256759722093
Current loss is: 0.06051962369139803
Current loss is: 0.060517567808465775
Current loss is: 0.06051702475103721
And now we plot the results,
error_ = u_predict .- u_realp1 = plot(xs,ys,u_real,linetype=:contourf,label = "analytic")p2 = plot(xs,ys,u_predict,linetype=:contourf,label = "predict")p3 = plot(xs,ys,error_,linetype=:contourf,label = "error")plot(p1,p2,p3)Integro Differential Equations
Let us consider an integro differential equation of the form:
formula not implementedwith boundary condition
formula not implemented t i(..)Di = Differential(t)Ii = Integral(t in DomainSets.ClosedInterval(0, t))eq = Di(i(t)) + 2*i(t) + 5*Ii(i(t)) ~ 1bcs = [i(0.) ~ 0.0]domains = [t ∈ Interval(0.0,2.0)]chain = Chain(Dense(1,15,Flux.σ),Dense(15,1))initθ = Float64.(DiffEqFlux.initial_params(chain))strategy_ = NeuralPDE.GridTraining(0.1)discretization = NeuralPDE.PhysicsInformedNN(chain, strategy_; init_params = nothing, phi = nothing, derivative = nothing, ) pde_system = PDESystem(eq,bcs,domains,[t],[i])prob = NeuralPDE.discretize(pde_system,discretization)res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=100)ts = [infimum(d.domain):0.01:supremum(d.domain) for d in domains][1]phi = discretization.phi┌ Warning: : no method matching get_unit for arguments (SciMLBase.NullParameters,).
└ @ ModelingToolkit /root/.julia/packages/ModelingToolkit/NRQ0J/src/systems/validation.jl:135
Current loss is: 2.324907204813763
Current loss is: 0.2528306155397627
Current loss is: 0.10577292220430932
Current loss is: 0.1051209376888921
Current loss is: 0.10158463217580985
Current loss is: 0.10042947006047223
Current loss is: 0.10000625290547008
Current loss is: 0.09907242390843148
Current loss is: 0.09832634598962825
Current loss is: 0.09683046741282958
Current loss is: 0.09526718978038413
Current loss is: 0.09241597622188294
Current loss is: 0.08803489569707795
Current loss is: 0.08379783318225129
Current loss is: 0.06388851179715133
Current loss is: 0.05678761896524551
Current loss is: 0.046856658339075694
Current loss is: 0.04234094624467016
Current loss is: 0.03887224280102542
Current loss is: 0.028522737873335965
Current loss is: 0.025143387883546617
Current loss is: 0.017849971934414576
Current loss is: 0.014076026067535875
Current loss is: 0.011157186145545317
Current loss is: 0.01073143632290276
Current loss is: 0.010194291463644515
Current loss is: 0.005177107483967417
Current loss is: 0.003253700397529548
Current loss is: 0.002611012449554509
Current loss is: 0.0025075887571848125
Current loss is: 0.002486161879581607
Current loss is: 0.00240414013510729
Current loss is: 0.001093106062313951
Current loss is: 0.00046286024097252985
Current loss is: 0.0002443257208143334
Current loss is: 0.00023733890816588608
Current loss is: 0.00023438208234169305
Current loss is: 0.00023384834700560694
Current loss is: 0.0002142997092826235
Current loss is: 0.0001792305101443119
Current loss is: 0.00012917645383039034
Current loss is: 0.00011184702248832961
Current loss is: 8.554851939880328e-5
Current loss is: 7.041827671936353e-5
Current loss is: 5.456540404768752e-5
Current loss is: 5.150984965501296e-5
Current loss is: 4.2556638278957334e-5
Current loss is: 3.802268097022839e-5
Current loss is: 3.335839647337833e-5
Current loss is: 3.29465397186763e-5
Current loss is: 3.162192757752035e-5
Current loss is: 3.134403614619635e-5
Current loss is: 3.1175545301838595e-5
Current loss is: 3.10376209779311e-5
Current loss is: 3.10309970402333e-5
Current loss is: 3.1019100954637576e-5
Current loss is: 3.099260336033678e-5
Current loss is: 3.0267188389514878e-5
Current loss is: 2.888275070492811e-5
Current loss is: 2.6540553288786213e-5
Current loss is: 2.46865366193367e-5
Current loss is: 2.1366177330336384e-5
Current loss is: 1.8255085755247105e-5
Current loss is: 1.4814527160884661e-5
Current loss is: 1.3626763840526283e-5
Current loss is: 1.3111792400868351e-5
Current loss is: 1.298075185035345e-5
Current loss is: 1.2963393639228268e-5
Current loss is: 1.2911586186937639e-5
Current loss is: 1.290365178851753e-5
Current loss is: 1.285446189292501e-5
Current loss is: 1.2833378972029226e-5
Current loss is: 1.259284862654859e-5
Current loss is: 1.2545785258522423e-5
Current loss is: 1.237231816632e-5
Current loss is: 1.2371623345662299e-5
Current loss is: 1.2328711451315788e-5
Current loss is: 1.2290063637453548e-5
Current loss is: 1.2180444163751699e-5
Current loss is: 1.2180195541066695e-5
Current loss is: 1.217004287650591e-5
Current loss is: 1.217004287650591e-5
Current loss is: 1.217004287650591e-5
Plotting the analytical solution and the predicted solution we get,
analytic_sol_func(t) = 1/2*(exp(-t))*(sin(2*t))u_real = [analytic_sol_func(t) for t in ts]u_predict = [first(phi([t],res.minimizer)) for t in ts]plot(ts,u_real, label = "Analytical Solution")plot!(ts,u_predict, label = "Predicted Solution")