Matrix PDEs

[seadra]

Any chance this NeuralPDE.jl can be eventually used with (complex) matrix PDEs?

Currently,

eq  = Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)*[0 1; 0 1]

gives an error

MethodError: no method matching ~(::Operation, ::Array{Operation,2})

Closest candidates are:

~(::Expression, !Matched::Expression) at /home/user/.julia/packages/ModelingToolkit/jYRzA/src/equations.jl:38

~(::Expression, !Matched::Number) at /home/user/.julia/packages/ModelingToolkit/jYRzA/src/equations.jl:39

~(::Number, !Matched::Expression) at /home/user/.julia/packages/ModelingToolkit/jYRzA/src/equations.jl:40

Stacktrace:

[1] top-level scope at In[17]:7

I have two questions:

• I assume this is a limitation ModelingToolkit, but maybe there an alternative way to define the PDEs with NeuralPDE.jl as a workaround?
• After a workaround or a proper fix in ModelingToolkit, will NeuralPDE.jl still require additional work to support a complex matrix PDE?

[ChrisRackauckas]

What do you want that equation to say? I don't understand the equation.

[seadra]

Sorry about the lack of details in the example. u is supposed to be a 4x4 matrix in the equation, although I couldn't find a way to specify it.

In practice, though, the kind of PDE I'm after is a complex vector equation (like Schrodinger equation with spin) where the Hamiltonian is a position-dependent matrix.

[ChrisRackauckas]

@variables u[1:4,1:4] and

A = sin(pi*x)*sin(pi*y)*[0 1; 0 1]
eq  = @. Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -A

?

[seadra]

Thanks for the response!
I just tried it, but it's giving an error:

MethodError: objects of type Array{Operation,2} are not callable
Use square brackets [] for indexing an Array.

Stacktrace:
 [1] _broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
 [2] _broadcast_getindex at ./broadcast.jl:621 [inlined]
 [3] _getindex at ./broadcast.jl:645 [inlined]
 [4] _broadcast_getindex at ./broadcast.jl:620 [inlined]
 [5] _getindex at ./broadcast.jl:644 [inlined]
 [6] _broadcast_getindex at ./broadcast.jl:620 [inlined]
 [7] _getindex at ./broadcast.jl:644 [inlined]
 [8] _broadcast_getindex at ./broadcast.jl:620 [inlined]
 [9] getindex at ./broadcast.jl:575 [inlined]
 [10] copy at ./broadcast.jl:876 [inlined]
 [11] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(~),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(+),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,Base.var"#62#63"{Differential,Differential},Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,Array{Operation,2},Tuple{Operation,Operation,Operation}}}},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,Base.var"#62#63"{Differential,Differential},Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,Array{Operation,2},Tuple{Operation,Operation,Operation}}}}}},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(-),Tuple{Array{Operation,2}}}}}) at ./broadcast.jl:837
 [12] top-level scope at In[9]:9

[ChrisRackauckas]

I don't think support for systems has been tagged yet, but @KirillZubov you might want to test on this one.

[KirillZubov]

Hi @seadra, the matrix PDEs form isn't supported for PINNs yet. But it will be implemented.

[ChrisRackauckas]

A matrix is just a system when you vec it, so once systems of PDEs are done this should be done.

[KirillZubov]

@seadra @ChrisRackauckas I'm sorry,
Actually, matrix PDE is already supported.
Here's an example that runs:

@parameters x y θ
@variables u[1:2,1:2](..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# matrix PDE
eqs  = @. [(Dxx(u_(x,y,θ)) + Dyy(u_(x,y,θ))) for u_ in u] ~ -sin(pi*x)*sin(pi*y)*[0 1; 0 1]

# Initial and boundary conditions
bcs = [[u[1](x,0,θ) ~ x u[2](x,0,θ) ~ 2;
        u[3](x,0,θ) ~ 3 u[4](x,0,θ) ~ 4]]
# or
# bcs = [@. [u_(x,0,θ) for u_ in u]  ~ [x 1; 3 4]]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0), y ∈ IntervalDomain(0.0,1.0)]

# Discretization
dx = 0.1
discretization = NeuralPDE.PhysicsInformedNN(dx)

# Neural network and optimizer
opt = Flux.ADAM(0.1)
chain = FastChain(FastDense(2,8,Flux.σ),FastDense(8,4))

pde_system = PDESystem(eqs,bcs,domains,[x,y],u)
prob = NeuralPDE.discretize(pde_system,discretization)
alg = NeuralPDE.NNDE(chain,opt,autodiff=false)
phi,res = solve(prob,alg,verbose=true, maxiters=50)

I'll double-check this task again when I move to generalize the parser.

[asdasdqwdqwfsdf]

@ChrisRackauckas can we test/run the "NeuralPDE.jl" examples in Binder Julia?

[ChrisRackauckas]

That should be possible.