Help me : PDE in Mathematica

Algebra, Applied Mathematics, Calculus, Differential Equations, Linear Algebra, Discrete Mathematics, Geometry, Number Theory, Complex Analysis, Differential Geometry, Dynamical Systems, Statistics, etc.
Forum Rules
By using the Wolfram Faculty Program Forum, you agree not to post any abusive, obscene, vulgar, slanderous, hateful, threatening, or sexually oriented material. Wolfram Faculty Program Forum administrators have the right to remove, edit, move or close any topic at any time should we see fit.

Personal Information: Posts in this forum may be viewed by non-members; however, the forum prohibits non-members from viewing your profile. Although your email address is hidden from both non-members and members, your account is initially configured to allow members to contact you via email through the forum. If you wish to hide your profile, or prohibit others from contacting you directly, you may change these settings by updating your profile through the User Control Panel.

Attachments: Attachments are not currently enabled on this forum. To share a file with others on this site, simply upload your file to the online storage service of your choice and include a link to the file within your post. If your school does not offer an online file storage and sharing service, the following sites provide free basic online file storage and sharing: Mozy, FilesAnywhere, Adrive, and KeepandShare.

Help me : PDE in Mathematica

Postby Dao_Trinh » Sun Jun 20, 2010 10:02 pm

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 !
dao.trinh@upmc.fr
User avatar
Dao_Trinh
 
Posts: 9
Joined: Thu May 06, 2010 3:22 pm
Location: Paris
Organization: Universite Pierre Marie Curie
Department: Chemistry

Re: Help me : PDE in Mathematica

Postby Kathy_Bautista » Thu Jun 24, 2010 12:05 am

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]
Katherine Bautista
Senior Academic Program Manager
Wolfram Research, Inc.
http://www.wolfram.com
User avatar
Kathy_Bautista
Site Admin
 
Posts: 182
Joined: Fri Jul 31, 2009 6:24 pm
Location: Mesa, Arizona
Organization: Wolfram Research, Inc.
Department: Academic Initiatives

Re: Help me : PDE in Mathematica

Postby Dao_Trinh » Thu Jun 24, 2010 11:33 am

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*
dao.trinh@upmc.fr
User avatar
Dao_Trinh
 
Posts: 9
Joined: Thu May 06, 2010 3:22 pm
Location: Paris
Organization: Universite Pierre Marie Curie
Department: Chemistry

Re: Help me : PDE in Mathematica

Postby Kathy_Bautista » Fri Jul 02, 2010 12:06 am

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
Katherine Bautista
Senior Academic Program Manager
Wolfram Research, Inc.
http://www.wolfram.com
User avatar
Kathy_Bautista
Site Admin
 
Posts: 182
Joined: Fri Jul 31, 2009 6:24 pm
Location: Mesa, Arizona
Organization: Wolfram Research, Inc.
Department: Academic Initiatives

Re: Help me : PDE in Mathematica

Postby Dao_Trinh » Sun Jul 04, 2010 2:36 pm

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
dao.trinh@upmc.fr
User avatar
Dao_Trinh
 
Posts: 9
Joined: Thu May 06, 2010 3:22 pm
Location: Paris
Organization: Universite Pierre Marie Curie
Department: Chemistry


Return to Mathematics and Statistics (Higher Education)

Who is online

Users browsing this forum: No registered users and 3 guests

cron