Help me : PDE in Mathematica

[Dao_Trinh]

Dear all,

I have a PDE in cylindrical coordinate but i cant imply the boundary condition :
In the text book, i have this equation :

eqn = D[u[r,z],{z,2}]+D[u[r,z],{r,2}+D[u[r,z],{r,1}]*1/r == 0

b1 = ( D[u[r,z],{z,1}]/.z->0 ) ==0
b2 = ( D[u[r,z],{r,1}]/.r->10^-5 ) ==0
b3= u[r,2]==1
b4 =u[2,z] ==1

NDSolve[{eqn,b1,b2,b3,b4},u,{r,0,2},{z,0,2}] ==> Error : u[2,z]==1 is not specified on single edge

and i dont use b4 :
NDSolve[{eqn,b1,b2,b3},u,{r,0,2},{z,0,2}] ===> Error : Number of constraint (1) is not equal total diff (2).

I cant find the solution ....

==================================================================================

In fact, i want to find the solution of this problem, using NDSolve. The problem is stated in this demonstration :
http://demonstrations.wolfram.com/Scann ... iskElectr/


Thanks for your help !

[Kathy_Bautista]

I have spoken to our technical support engineers about this issue, and they commented that this is basically related to how initial/boundary conditions are to be specified for NDSolve. I have attached their comments below.

-Kathy

-----------------
Certain limitations exist in terms of how these can be specified. I have tried to give some workarounds to get a reasonable solution, which you can see below.

If the problem (system) has singularities, then PrecisionGoal may have to be used to get better results.

Given equation and boundary conditions

IN:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] + D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 10^-5) == 0; b3 = u[r, 2] == 1; b4 =
u[2, z] == 1;
NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 0, 2}, {z, 0, 2}]

ERROR:
NDSolve::bcedge: Boundary condition u[2,z]==1 is not specified on a single edge of the boundary of the computational domain. >>

The above error indicates that there is inconsistency between the specified boundary conditions and the domain of NDSolve. Note that the domain of NDSolve is given as {r,0,2} and {z,0,2}.

Try changing the b.c to fit the domain values for r.

IN:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 10^-5) == 0; b3 =
u[r, 2] == 1; b4 = u[10^-5, z] == 1;
NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 10^-5, 2}, {z, 0, 2}]

ERROR:
NDSolve::eerr: Warning: Scaled local spatial error estimate of 717.6034425350999` at r = 2.` in the direction of independent variable z is much greater than prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or you may want to specify a smaller grid spacing using the MaxStepSize or MinPoints method options. >>

The above warning message indicates that there could be some singularity in the solution, and the error can be high.

Try using PrecisionGoal

IN:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 10^-5) == 0; b3 =
u[r, 2] == 1; b4 = u[10^-5, z] == 1;
NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 10^-5, 2}, {z, 0, 2},
PrecisionGoal -> 2]

The solution looks reasonable.

Notice that NDSolve must have "consistent" boundary and initial conditions. Thats why the following fails. Note the error message.

IN:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 10^-5) == 0; b3 =
u[r, 2] == 1; b4 = u[2, z] == 1;
NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 10^-5, 2}, {z, 0, 2},
PrecisionGoal -> 2]

ERROR:
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. >>

The following also works. (Changing the boundary conditions for r.)

IN:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 2) == 0; b3 = u[r, 2] == 1; b4 =
u[2, z] == 1;
NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 10^-5, 2}, {z, 0, 2},
PrecisionGoal -> 2]

[Dao_Trinh]

Dear Kathy Bautista,

I feel really happy with the solution. Thank you very much Kathy Bautista and the support of technical team.

That is the first step i reached. I know how to apply the homogeneous boundary in my problem.

The second step is how to imply the mixed boundary in my problem (as stated in the demonstration http://demonstrations.wolfram.com/Scann ... iskElectr/.

This mixed boundary is :

For z = 0 :
When 0 < r < 1 : we used the Dirichlet boundary : u[ r,0] == 0
when 1 < r < 2 : we used the Neumann boundary : du/dz ==0

So, how to imply this mixed boundary in Mathematica for 0 < r < 2. I searched in Wolfram library but i found nothing about the mixed boundary or Robin boundary.

I tried to separate the problem in 2 domains, for u[r,z] r < 1 and u[r,z] r > 1. But it seems that the boundary isn't the same in 2 domains.

My problem is a difficult problem. But i found that in the text book and article *, they can solve my problem by Finite Difference, and i want to try to solve it with NDSolve. Also, i have the analytical solution (stated in the demonstration).

Best regards and thanks for your precious help

Dao Trinh


* The article which stated my problem with boundary condition:
Journal of Electroanalytical Chemistry
Volume 456, Issues 1-2, 30 September 1998, Pages 1-12
D. J. Gavaghan*

[Kathy_Bautista]

My colleague Aravind Hanasoge sent me notes on one way to approach this problem. I have attached his comments below...

-Kathy

-----------------

Basically, the equations are solved for two domains of r. In the first domain, 0<r<1, and one of the boundary conditions used is u[r,0]==0.

In the second domain, 1<r<2, one of the boundary conditions used is du/dz|z=0 == 0

The other boundary conditions are left unchanged.

We get two interpolation functions as two solutions. Then we can construct a Piecewise function based on these two solutions to get a complete solution.

You may want to change the other boundary conditions to suit your needs.

Input:
eqn1 = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r == 0; b11 = (u[r, z] /. z -> 0) ==
0; b21 = (D[u[r, z], {r, 1}] /. r -> 10^-5) == 0; b31 =
u[r, 2] == 1; b41 = u[10^-5, z] == 1;
u1 = NDSolve[{eqn1, b11, b21, b31, b41}, u, {r, 10^-5, 1}, {z, 0, 2},
PrecisionGoal -> 2][[1, 1, 2]]

Output:
NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent. >>
InterpolatingFunction[{{0.00001, 1.},{0., 2.}},<>]

Input:
eqn = D[u[r, z], {z, 2}] + D[u[r, z], {r, 2}] +
D[u[r, z], {r, 1}]*1/r ==
0; b1 = (D[u[r, z], {z, 1}] /. z -> 0) ==
0; b2 = (D[u[r, z], {r, 1}] /. r -> 1) == 0; b3 = u[r, 2] == 1; b4 =
u[1, z] == 1;
u2 = NDSolve[{eqn, b1, b2, b3, b4}, u, {r, 1, 2}, {z, 0, 2},
PrecisionGoal -> 2][[1, 1, 2]]

Output:
InterpolatingFunction[{{1., 2.}, {0., 2.}},<>]

Input:
usol = Piecewise[{{u1, 0 < r < 1}, {u2, 1 <= r < 2}}]

Output:
InterpolatingFunction[{{0.00001, 1.}, {0., 2.}},<>] 0 < r < 1
InterpolatingFunction[{{1., 2.}, {0., 2.}},<>] 1 <= r < 2
0 True

[Dao_Trinh]

Dear Kathy_Bautista and Aravind Hanasoge,

Thank you very much for your help in my problem.

Now, i know how to combine the boundary in 2 domains. But the solution (of NDSolve) is not the same with the analytical solution, i don't know why.

Anyway, i appreciate so much your help, thank you again.

Best regards,
Dao Trinh