EDIT: actualy the nonlinear partial differential equations for interacting density distributions, including boundary conditions, should be given as

$$ \frac{\partial\phi}{\partial t} = D \frac{\partial^{2}\phi}{\partial x^{2}}, \ \ \quad \qquad \frac{\partial\phi_{1}}{\partial t} = D \frac{\partial^{2}\phi_1}{\partial x^{2}} +\frac{\partial^{2}}{\partial x^{2}}(\frac{\phi_{1}}{\phi})\\ \frac{\partial\phi(t,-1)}{\partial x}=\frac{\partial\phi(t,1)}{\partial x}=0, \ \ \text{and} \qquad \frac{\partial\phi_1(t,-1)}{\partial x}=\frac{\partial\phi_1(t,1)}{\partial x}=0. $$

I am trying to solve the following nonlinear partial differential equations,

$$ \frac{\partial\phi}{\partial t} = D \frac{\partial^{2}\phi}{\partial x^{2}}, \ \ \text{and} \qquad \frac{\partial\phi_{1}}{\partial t} = D \frac{\partial^{2}\phi_1}{\partial x^{2}} -\frac{\partial^{2}}{\partial x^{2}}(\frac{\phi_{1}}{\phi}). $$

Where both $\phi(t,x)$ and $\phi_1(t,x)$ are $ \mathbb{R}^2 \rightarrow \mathbb{R}_{>0}$ functions of $t,\ x$. With boundary conditions,

$$ \frac{\partial\phi(t,-1)}{\partial x}=\frac{\partial\phi(t,1)}{\partial x}=0, \ \ \quad \qquad \frac{\partial\phi_1(t,-1)}{\partial x}=\frac{\partial\phi_1(t,1)}{\partial x}=0, \\ \frac{\partial\phi(0,x)}{\partial t}=0,\ \ \text{and} \qquad \frac{\partial\phi_1(0,x)}{\partial t}=0. $$

The problem above describes two density distributions $\phi_1$, $\phi_2$ and their sum $\phi$, that diffuse over time along the $x$-axis between $x=-1$ and $x=1$. Furthermore, the assumsion is made that while the sum of the distributions diffuses normaly, $\phi_1$ and $\phi_2$ will interact in such a way that the flux of the first distribution, from the areas with the high ratio $\frac{\phi_1}{\phi}$ to the areas with lower values of this ratio, is given by $\frac{\partial}{\partial x} (\frac{\phi_1}{\phi})$.

For my problem, the initial conditions are given by two normal distributions with their mean values between $x=-1$ and $x=1$, idealy I would like to find a way of numericaly solve the pde's for any such initial conditions.

My mathematica code for solving just $\phi$ works like a charm, with no errors whatsoever.

pde = {D[ϕ[t, x], t] == D[ϕ[t, x], x, x]};

ic = {ϕ[0, x] == PDF[NormalDistribution[-.3, .1], x] + PDF[NormalDistribution[.1, .2], x]}; 

bc = {
Derivative[1, 0][ϕ][0, x] == 0,
Derivative[0, 1][ϕ][t, -1] == 0,
Derivative[0, 1][ϕ][t, 1] == 0};

sol = NDSolve[Join[pde, ic, bc], ϕ, {t, 0, 1}, {x, -1, 1}]

Plot3D[{ϕ[t, x] /. sol[[1]]}, {t, 0, .1}, {x, -1, 1}, PlotRange -> All]

Using this anwer as a numerical solution and solving only for $\phi_1$ works marginally, since (I think)

ϕ[t, x] /. sol[[1]]

isn't always well defined, such that the quotient $\frac{\phi_1}{\phi}$ becomes very large and messes up the solotion.

pde2 = {
D[Subscript[ϕ, 1][t, x], t] == D[Subscript[ϕ, 1][t, x], x, x] - D[Subscript[ϕ, 1][t, x]/(ϕ[t, x] /. sol[[1]]), x, x]};

ic2 = {Subscript[ϕ, 1][0, x] == PDF[NormalDistribution[.1, .2], x]};

bc2 = {
Derivative[1, 0][Subscript[ϕ, 1]][0, x] == 0,
Derivative[0, 1][Subscript[ϕ, 1]][t, -1] == 0,
Derivative[0, 1][Subscript[ϕ, 1]][t, 1] == 0};

sol2 = NDSolve[Join[pde2, ic2, bc2], Subscript[ϕ, 1], {t, 0, .1}, {x, -1, 1}];

Plot3D[{Subscript[ϕ, 1][t, x] /. sol2}, {t, 0, .1}, {x, -1, 1}]

But when I try to simultaneously solve for $\phi$ and $\phi_1$ errors appear.

EDIT: the code below is edited and works now. The error originated in the improper partial differential equations due to a faulty assumption and a redundent boundary condition.

pde3 = {
D[ϕ[t, x], t] == D[ϕ[t, x], x, x],
D[Subscript[ϕ, 1][t, x], t] == D[Subscript[ϕ, 1][t, x], x, x] + D[Subscript[ϕ, 1][t, x]/ϕ[t, x], x, x]};

ic3 = {
ϕ[0, x] == PDF[NormalDistribution[-.3, .1], x] + PDF[NormalDistribution[.1, .2], x],
Subscript[ϕ, 1][0, x] == PDF[NormalDistribution[.1, .2], x]};

bc3 = {
Derivative[0, 1][ϕ][t, -1] == 0,
Derivative[0, 1][ϕ][t, 1] == 0,

Derivative[0, 1][Subscript[ϕ, 1]][t, -1] == 0,
Derivative[0, 1][Subscript[ϕ, 1]][t, 1] == 0};

sol3 = NDSolveValue[Join[pde3, ic3, bc3], {ϕ, Subscript[ϕ, 1]}, {t,0, .1}, {x, -1, 1}]

So far I believe the error's originate from the inproper method used by NDsolve, but I am not sufficiently skilled in this area to solve this on my own. Then again I could be something else entirly. Any help would be greatly appreciated.


closed as off-topic by xzczd, bbgodfrey, Sjoerd C. de Vries, Michael E2, Dr. belisarius May 6 '15 at 18:09

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – xzczd, bbgodfrey, Sjoerd C. de Vries, Michael E2, Dr. belisarius
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ v10 manages to handle bc[[1]] and bc2[[1]]? I think they're redundant, actually they cause failures of NDSolve in v9. $\endgroup$ – xzczd May 5 '15 at 12:24
  • $\begingroup$ Indeed I forgot to mention I am on v10. And yes, v10 manages to handle both bc[[1]] and bc[[2]]. Omitting them results in, for atleast $\phi$, the same solution. I added the last boundary conditions when I was looking for movable endpoints on $x=-1$ and $x=1$ (the source was for the wave equation), so I wasn't sure if they where needed. But when omitted NDSolve gives a warning that het boundary and inital conditions are inconsistent, so I kept them. But now seeing they are realy inconsistent at $x=-1$ and $x=1$, I am inclined to agree. $\endgroup$ – user19218 May 5 '15 at 12:44
  • $\begingroup$ @xzczd Omitting bc[[1]] and bc[[2]], then sequentially evaluating sol1 and sol2 yields a result for sol1 almost instantaneously, but sol2 seems to keep running indefinitely and using alot of CPU. $\endgroup$ – user19218 May 5 '15 at 13:12
  • $\begingroup$ @xzczd I redid the math and concluded I made a mistake, the flux from the distribution should be negatively proportional to $\frac{\partial}{\partial x}\frac{\phi_1}{\phi}$. This change, and omitting the boundary conditions earlier mentioned, leads to a solution without any hickups. Thank you for the fresh pair of eyes, this helped alot! $\endgroup$ – user19218 May 5 '15 at 14:02