# How to solve PDAEs (incorporate mass matrix into PDE solver)

I am trying to solve a system of 1st order linear partial differential equations with a constraint on two out of three of my dependent variables (to enforce population conservation).

I am familiar with DAEs and the associated documentation, and I have built my PDE solver by discretizing my spatial derivative using the first method outlined here (Solving Large Stiff Equations · DifferentialEquations.jl), but I have been struggling to implement the mass matrix (or the implicit DAE approach) into my PDE solver.

My first attempt at a mass-matrix approach (with the dimensions num_dependent_vars x num_dependent_vars x num_spatial_points) raised a dimensions-related error in my code (mass matrix dimensions must match the length of my initial condition “vector”, which is num_spatial_points x num_dependent_vars).

My question is: how should I incorporate an algebraic constraint on my dependent variables into my semi-discretized PDE solver? I haven’t been able to find any documentation relevant to this specific problem. I apologize in advance if I have just been missing the obvious. I don’t have much experience creating MWEs but I have attempted to simplify my code into a mostly relevant MWE which produces an error similar (although, bafflingly, not identical) to the one I stated above.

Please let me know if there is anything I can do to clarify this question. I’m also not sure if this is quite the right place to ask it but I was not sure where else to do so. Any assistance or advice would be greatly appreciated!

“MWE”:

``````using DifferentialEquations, LinearAlgebra

xgrid = collect(-pi:0.2:pi)

limit(a,N) = a == N+1 ? 1 : a == 0 ? N : a

# start list of functions
function efield(t,Emax,sig,w0)

E_t = Emax * exp(-4 * log(2) * (t/sig)^2) * cos(t * w0) ;

return E_t
end

function pdesystem(dx,x,xgrid,t)

u = @view x[:,1] ; # soc
v = @view x[:,2] ;
w = @view x[:,3] ;

knum = length(xgrid)
T2 = 500
Eg = 1.5 / 27.21 ; # eV -> AU conversion
Emax = 0.0014;
sig = 2687;
w0 = 0.006;
Rabi = 9.5 * efield(t,Emax,sig,w0);
dy = xgrid - xgrid; # for evenly spaced 1D line of x points
for i = 1:knum
ip1, im1 = limit(i+1,knum), limit(i-1,knum)
dx[i,1] = ( E * (u[ip1] - u[im1])/2/dy ) - u[i] / T2  +  Eg[i] * v[i] ;
dx[i,2] = ( E * (v[ip1] - v[im1])/2/dy ) - Eg[i] * u[i]  -  2 * Rabi * w[i] - v[i] / T2;
dx[i,3] = ( E * (w[ip1] - w[im1])/2/dy ) + 2 * Rabi * v[i] ;
dx[i,4] = 1 - 2 * w[i] ;
end

end

function init_sbe(xgrid)

N = length(xgrid);
u = Array{Float64}(undef,(N,4));
u[i,1] = 0
u[i,2] = 0
u[i,3] = -1
u[i,4] = 0
end
return u
end

function massmx(xgrid)
N = length(xgrid);
M = Array{Float64}(undef,(N,4,4));
fill!(M,0)
M[i,1,1] = 1
M[i,2,2] = 1
M[i,3,3] = 1
M[i,4,4] = 0
end
return M
end

# problem set up
u0 = init_sbe(xgrid);
Nk = length(xgrid);
trange = (-12396.694214876032, 12396.694214876032) ;
M = massmx(xgrid);
f1 = ODEFunction(pdesystem,mass_matrix=M)
prob_sbe = ODEProblem(f1,u0,trange,xgrid) ;
soln1 = solve(prob_sbe,Rodas5(),reltol=1e-10);

``````

The mass matrix should be a matrix whose properties match the `vec(u)` form of the equation. so for `u` sized `Nx4`, you treat this as `vec(u)` of length `4N` and so it should be a `4N x 4N` matrix. Then for your case, the top left 3N x 3N matrix will be diagonal, and the other “quadrants” would be 0 matrices

Thank you so much for your response!! This fixes my MWE! I did not realize the ODEFunction would handle the necessary dimensions reshaping automatically, this is very convenient!

When I implement this fix in my actual code, I get ERROR: Singular Exception (1) or ERROR: Singular Exception (18010), so I just want to verify: this means there is some instability in the solver, correct? And thus likely an issue with my model/solver algorithm/tolerance settings rather than a problem with the way I have set up the DAE?

Is your DAE index 1? I.e., are all of your constraints functions of the algebraic variables?

In my actual code (not the mwe), I have several dependent variables, and my constraint on three of those dependent variables is (basically): w(x,t) - u(x,t) - v(x,t) = 0 (i.e., u(x,t)+ v(x,t) = w(x,t)) which I believe is of index 1. Apologies if I am incorrect, I am rather new to all this…

One of those dependent variables is not a differential variable, i.e. it has no equation where it shows up having a derivative in time?