Periodic boundary conditions with multiple variables

Question

I am trying to numerically solve the following first order coupled differential equations numerically, where i is an integer (can be set to zero), le = 1, lb = 1, c = 0.5, and Ly is taken to be 10. I am receiving the error that

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

It seems that no matter how I try to impose periodicity on the solution it returns the error. Does someone know how to impose periodic boundary conditions on this type of PDE?

i = 0; le = 1; lb = 1; c = 0.5; Ly = 10;

NDSolve[{(-I)*D[ψ1[x, y, z, t], t] + (c + 1)*((-I)*D[ψ1[x, y, z, t], z] - (t*ψ1[x, y, z, t])/le^2) + 
     ((-I)*D[ψ2[x, y, z, t], x] + (y/lb^2)*ψ2[x, y, z, t] - D[ψ2[x, y, z, t], y]) == 10^(-3), 
   (-I)*D[ψ2[x, y, z, t], t] + (c - 1)*((-I)*D[ψ2[x, y, z, t], z] - (t*ψ2[x, y, z, t])/le^2) + 
     (-I)*D[ψ1[x, y, z, t], x] + (y/lb^2)*ψ1[x, y, z, t] + D[ψ1[x, y, z, t], y] == 10^(-3), 
   PeriodicBoundaryCondition[ψ1[x, y, z, t], y == -((Ly*((100 - i)/100))/2), 
    Function[y, y + Ly*((100 - i)/100)]], PeriodicBoundaryCondition[ψ2[x, y, z, t], 
    y == -((Ly*((100 - i)/100))/2), Function[y, y + Ly*((100 - i)/100)]], 
   PeriodicBoundaryCondition[ψ1[x, y, z, t], x == -((2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100)))), 
    Function[x, x + (2*Pi*20)/(lb^2*(Ly*((100 - i)/100)))]], PeriodicBoundaryCondition[ψ2[x, y, z, t], 
    x == -((2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100)))), 
    Function[x, x + (2*Pi*20)/(lb^2*(Ly*((100 - i)/100)))]], PeriodicBoundaryCondition[ψ1[x, y, z, t], 
    t == (Ly*((100 - i)/100))/2, Function[t, t - Ly*((100 - i)/100)]], 
   PeriodicBoundaryCondition[ψ2[x, y, z, t], t == (Ly*((100 - i)/100))/2, 
    Function[t, t - Ly*((100 - i)/100)]], PeriodicBoundaryCondition[ψ1[x, y, z, t], 
    z == (2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100))), Function[z, z - (2*Pi*20)/(lb^2*(Ly*((100 - i)/100)))]], 
   PeriodicBoundaryCondition[ψ2[x, y, z, t], z == (2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100))), 
    Function[z, z - (2*Pi*20)/(lb^2*(Ly*((100 - i)/100)))]]}, {ψ1[x, y, z, t], ψ2[x, y, z, t]}, 
  {x, (-2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100))), (2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100)))}, 
  {y, ((-Ly)*((100 - i)/100))/2, (Ly*((100 - i)/100))/2}, {z, (-2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100))), 
   (2*Pi*20)/(lb^2*(2*Ly*((100 - i)/100)))}, {t, ((-Ly)*((100 - i)/100))/2, (Ly*((100 - i)/100))/2}, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}}]

Answer 1

You'd need to add initial conditions like

\[Psi]1[x, y, z, 0] == 0, \[Psi]2[x, y, z, 0] == 0

Then, the method option should be something like:

Method -> {"PDEDiscretization" -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.1}}}}

You could also consider not using the FEM and have periodic boundary conditions of the form

\[Psi]1[xLeft, y, z, t] == \[Psi]1[xRight, y, z, t]