# Symbolic ODE/PDE in Vector form with Symbolics/ModelingToolkit

**URL:** https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793
**Category:** Modelling & Simulations
**Tags:** pde, ode, modelingtoolkit, symbolics, modelling
**Created:** [September 12, 2023, 5:33pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793 "2023-09-12T17:33:48Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![Bizzi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bizzi/32/51484_2.png) [@Bizzi](https://discourse.julialang.org/u/Bizzi)
#### Post date: [September 12, 2023, 5:33pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/1 "2023-09-12T17:33:48Z")

</div>

I’ve been having a hard time trying to use _vector-valued equations_ (as opposed to vectors of equations) in Symbolics/ModelingToolkit. I remember reading somewhere in either docs that his feature was somewhat experimental, yet I’ve been unable to find this piece of information again, as well as any tutorial or guidance on how to do this.

Should I just leave this alone? My work involves large ODEs in vector/matrix form and my solver (a PINN) gives me solutions also in vector form; It just feels wrong to have to split both the equation and the solution into their individual lines and link them one by one.

An example of the libraries not working how one (I?) would expect:

```julia
using ModelingToolkit

@variables t, u(..), v(..), U(..)[1:2]

#Simple Harmonic Oscillator
Dt = Differential(t) 
eqs = [Dt(u(t)) ~ v(t), Dt(v(t)) ~ -u(t)] #Works
vec_eqs = Dt(U(t)) ~ [0 1;-1 0]*U(t) #ERROR: axes of Differential(t)((U(t))[1:2]) not known

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 12, 2023, 5:52pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/2 "2023-09-12T17:52:53Z")

</div>

It follows Julia’s vector semantics, so you need to broadcast that:

```julia
@variables t, u(..), v(..), U(..)[1:2]
Dt = Differential(t) 
vec_eqs = Dt.(U(t)) .~ [0 1;-1 0]*U(t)

```

That could probably be simplified though.

---

<div class="post-metadata">

### Author: ![Bizzi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bizzi/32/51484_2.png) [@Bizzi](https://discourse.julialang.org/u/Bizzi)
#### Post date: [September 12, 2023, 7:48pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/3 "2023-09-12T19:48:27Z")

</div>

Thanks Chris. I had tried broadcasting, but it seems to mess up the types. In particular, `vec_eqs` is not an `Equation`, so it cannot be used with DifferentialEquations/NeuralPDE. Or am I missing something?

```julia
using ModelingToolkit

@parameters σ ρ
@variables t x(t) y(t) X(t)[1:2]
D = Differential(t)

eqs = [D(x) ~ σ*y,
       D(y) ~ ρ*x]
vec_eqs = D.(X) .~ [0 σ;ρ 0]*X

typeof(eqs)#Vector{Equation}
typeof(vec_eqs)#Symbolics.Arr{Any,1}

@named sys = ODESystem(eqs,t,[x,y],[σ,ρ],tspan=(0, 1000.0)) #Works
@named vec_sys = ODESystem(vec_eqs,t,[X],[σ,ρ],tspan=(0, 1000.0))#ERROR

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 13, 2023, 3:12am UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/4 "2023-09-13T03:12:58Z")

</div>

@xtalax were you working on vector derivatives in NeuralPDE?

---

<div class="post-metadata">

### Author: ![Bizzi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bizzi/32/51484_2.png) [@Bizzi](https://discourse.julialang.org/u/Bizzi)
#### Post date: [September 13, 2023, 12:36pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/5 "2023-09-13T12:36:42Z")

</div>

Please correct me if I’m saying something very silly, but it seems like NeuralPDE itself doesn’t really care if it’s dealing with vectors or scalars? As long as the equations themselves are understood by Symbolics things seem to work.

For example, take this (trivial) system of two conservation laws with a single PINN:

```julia
#Vector_Neural_PDE
using NeuralPDE, Lux, Optimization, OptimizationOptimJL
import ModelingToolkit: Interval

@parameters x t
@variables u(..) #Declare u as if it were a scalar
Dx = Differential(x)
Dt = Differential(t)

#Set PDE normally
c = 1
eq = Dt(u(x,t)) ~ c*Dx(u(x,t))

#Give a vector disguised as a scalar for boundary conditions
bcs = [u(x,0) .~ exp(-x^2)] 

# Space and time domains
domains = [x ∈ Interval(-5,5),
           t ∈ Interval(0.0,10.0)]

# PINN with TWO outputs
dim_in = 2
dim_out = 2
chain = Lux.Chain(Dense(dim_in,12,Lux.tanh),Dense(12,12,Lux.tanh),Dense(12,12,Lux.tanh),Dense(12,dim_out))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(chain,GridTraining(dx))
@named pde_system = PDESystem(eq,bcs,domains,[x,t],[u(x, t)])
prob = discretize(pde_system,discretization)

#Optimizer
opt = OptimizationOptimJL.BFGS()

res = Optimization.solve(prob, opt, maxtime=20)

```

All of this seems to work just fine, with both outputs of our PINN being trained accordingly:

```julia
using Plots

phi = discretization.phi

xs,ts = [infimum(d.domain):dx/10:supremum(d.domain) for d in domains]

U(x,t) = phi([x,t],res.u)[1]
V(x,t) = phi([x,t],res.u)[2]

plot(xs,U.(xs,0))
plot!(xs,U.(xs,1))
plot!(xs,U.(xs,2))

plot(xs,V.(xs,0))
plot!(xs,V.(xs,1))
plot!(xs,V.(xs,2))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/8/f80b8734716f843e86a998c7cd413e6556b772ee.png)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/1/b1b34d6c272ba42cd8a9a91b938d94d7ddd12e1f.png)

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 14, 2023, 11:26am UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/6 "2023-09-14T11:26:20Z")

</div>

> [@Bizzi](#):
>
> Please correct me if I’m saying something very silly, but it seems like NeuralPDE itself doesn’t really care if it’s dealing with vectors or scalars? As long as the equations themselves are understood by Symbolics things seem to work.

Yes, that works. But the training process could specialize on vector calculus to improve its operations.

---

<div class="post-metadata">

### Author: ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)
#### Post date: [September 14, 2023, 1:52pm UTC](https://discourse.julialang.org/t/symbolic-ode-pde-in-vector-form-with-symbolics-modelingtoolkit/103793/7 "2023-09-14T13:52:06Z")

</div>

Yes, see here [Add ArrayDifferentialOperators for Vector calculus by xtalax · Pull Request #942 · JuliaSymbolics/Symbolics.jl · GitHub](https://github.com/JuliaSymbolics/Symbolics.jl/pull/942)
