# How to solve Coupled a Parabolic and Elliptic PDE in NDSolve?

I want to solve a mixed PDE Parabolic-Elliptic system in 3-dimension (rectangular coordinate), as shown below:

The respective code version with parameters value, boundary and initial conditions is,

    L = 1000;(*length of cube*)
pts = 200;
T = 400;(*Time integration*)
ϵ = 50 λ;
s = 2.75;
δ = 2.76;
γ = 2.75;
τ = 3.65;
χ = 1.10;
A = 1.6438;
μ = 0.2;
θ = 1;
λ = 50;

(*system of nonlinear PDE*)

pde = {0 == ϵ Laplacian[σ[t,x, y, z], {x, y, z}] + s - δ ϕ[t, x, y, z] - γ σ[t,x, y, z], D[ϕ[t, x, y, z],t] == λ Laplacian[ϕ[t, x, y, z], {x, y, z}] + 64/τ (1 - ϕ[t, x, y, z]) (ϕ[t, x, y, z] - 1/2) + χ σ[t,x, y, z] ϕ[t, x, y, z] - Aϕ[t, x, y, z] - 3 μ θ^2 (2 θ - 3) ϕ[t, x, y, z] (ϕ[t, x, y, z] - 1)};

(*Periodic boundary condition*)

bc = {ϕ[t, 0, y, z] == ϕ[t, L, y, z], ϕ[t, x, 0, z] == ϕ[t, x, L, z], ϕ[t, x, y, 0] == ϕ[t, x, y,L], σ[t, 0, y, z] == σ[t, L, y, z], σ[t, x, 0, z] == σ[t, x, L, z], σ[t, x, y, 0] == σ[t, x, y, L]};

(*initial condition*)

ic = {ϕ[0, x, y, z] == If[(x - 500)^2 + (y - 500)^2 + (z - 500)^2 <= (25)^2, 1, 0]};

eqns = Flatten@{pde, bc, ic};

(*Integration*)

sol = NDSolve[eqns, {ϕ, σ}, {t, 0, T}, {x, 0, L}, {y, 0, L}, {z, 0, L}, Method -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> pts, "MaxPoints" -> pts}}];


but something is not working well.

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.

Can anybody help me?

• First equation is unclear formulated. Do you suppose that $\frac{\partial \sigma}{\partial t}=0$? Nov 1, 2021 at 3:17
• @AlexTrounev , it is assumed the $\sigma$ to be in a stationary state in the time scale.
– SAC
Nov 1, 2021 at 3:51
• Ok, I understand that you try to solve pde in a cube. But why do you set L=1000 with pts=200? Is it real problem or you just test NDSolve? Nov 1, 2021 at 4:26
• @AlexTrounev, it is a real problem of mathematical modelling.
– SAC
Nov 1, 2021 at 11:36
• It looks like the Cahn-Hilliard model discussed on mathematica.stackexchange.com/questions/202446/… Nov 1, 2021 at 15:17

You can not solve mixed time dependent and stationary equations. You'd have to make the first equation time dependent. Something like this starts to time integrate, though I did not wait for it to finish:

pde = {D[σ[t, x, y, z], t] == ϵ Laplacian[σ[t, x, y, z], {x, y, z}] +
s - δ ϕ[t, x, y, z] - γ σ[t, x, y, z],
D[ϕ[t, x, y, z], t] == λ Laplacian[ϕ[t, x, y, z], {x, y, z}] +
64/τ (1 - ϕ[t, x, y, z]) (ϕ[t, x, y, z] -
1/2) + χ σ[t, x, y, z] ϕ[t, x, y, z] - A ϕ[t, x, y, z] -
3 μ θ^2 (2 θ - 3) ϕ[t, x, y, z] (ϕ[t, x, y, z] - 1)};

ic = {σ[0, x, y, z] == 0, ϕ[0, x, y, z] ==
If[(x - 500)^2 + (y - 500)^2 + (z - 500)^2 <= (25)^2, 1, 0]};

• @Nasser, thanks for the edit! Nov 2, 2021 at 10:57
• Thanks for the suggestion. I tried to integrate the equations as they are and it takes a long time. Does anyone have any idea how to optimize integration time? Could the Fast Fourier Transform (FFT) decrease simulation time?
– SAC
Nov 2, 2021 at 17:28