I'm trying to numerically solve a second-order partial differential equation in 4 dimensions. The PDE is a Fokker-Planck equation and the variables are the real and imaginary parts of the field of two coupled oscillators. I try to solve this with NDSolve using Dirichlet boundaries, but I get a strange error. I first create a 4-dimensional region in which the solution has to be found.

xmin = 0;
xmax = 13;
ymin = -10;
ymax = 10;
dbc = {DirichletCondition[P[x1, y1, x2, y2] == 0.0, True]};
Ω4d = 
 Parallelepiped[{xmin, ymin, xmin, 
   ymin}, {{xmax - xmin, 0, 0, 0}, {0, ymax - ymin, 0, 0}, {0, 0, 
    xmax - xmin, 0}, {0, 0, 0, ymax - ymin}}]

Then I create a very nasty PDE which I call pde $$\left(-0.02 \text{x1}^2 \text{y1}+1. \text{x1}^3+1. \text{x1} \text{y1}^2+1. \text{x1}+0.1 \text{x2}-0.02 \text{y1}^3-0.002 \text{y2}\right) P^{(0,1,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(-1. \text{x1}^2 \text{y1}-0.02 \text{x1}^3-0.02 \text{x1} \text{y1}^2-0.002 \text{x2}-1. \text{y1}^3-1. \text{y1}-0.1 \text{y2}\right) P^{(1,0,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(0.1 \text{x1}-0.02 \text{x2}^2 \text{y2}+1. \text{x2}^3+1. \text{x2} \text{y2}^2+1. \text{x2}+0.018 \text{y1}-0.02 \text{y2}^3-0.02 \text{y2}\right) P^{(0,0,0,1)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(0.018 \text{x1}-1. \text{x2}^2 \text{y2}-0.02 \text{x2}^3-0.02 \text{x2} \text{y2}^2-0.02 \text{x2}-0.1 \text{y1}-1. \text{y2}^3-1. \text{y2}\right) P^{(0,0,1,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+\left(-0.04 \text{x1}^2-0.04 \text{x2}^2-0.04 \text{y1}^2-0.04 \text{y2}^2-0.04\right) P(\text{x1},\text{y1},\text{x2},\text{y2})=0.02 P^{(0,0,0,2)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.02 P^{(0,0,2,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.002 P^{(0,2,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})+0.002 P^{(2,0,0,0)}(\text{x1},\text{y1},\text{x2},\text{y2})$$ I try to solve this with

NDSolve[{pde, dbc}, P, {x1, y1, x2, y2} ∈ Ω4d]

But I get the following error message

NDSolve::fem3: -- Message text not found -- (Parallelepiped[{0,-10,0,-10},{{13,0,0,0},{0,20,0,0},{0,0,13,0},{0,0,0,20}}])

If I solve with

NDSolve[{pde, dbc}, P, {x1, xmin, xmax}, {y1, ymin, ymax}, {x2, xmin,xmax}, {y2, ymin, ymax}]

I get

NDSolve::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.

Does anyone know what I'm doing wrong here?


1 Answer 1


The finite element method as implemented in Version 12.1 (and 12.2) can solve PDEs in space dimensions, 1, 2 and 3. So a 4D spatial PDE can not be solved with the FEM.

You can use the TensorProductGrid method to solve 4D PDEs but then you will have to use a rectangular domain.

  • $\begingroup$ Thanks! Do you know which method can solve a 4D PDE? $\endgroup$
    – Philip
    Dec 10, 2020 at 15:24
  • $\begingroup$ @Philip, see update. $\endgroup$
    – user21
    Dec 10, 2020 at 17:09
  • $\begingroup$ All numerical methods such as FD, RBF-FD and RBF-HFD can be used for solving a 4D PDE problem on TensorProductGrid. Especially their variants on non-uniform grids might be better. However the solvers may suffer of curse of dimensionality and thus a careful tuning of some input variables might be needed. $\endgroup$
    – Faz
    Jul 6, 2023 at 14:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.