5
$\begingroup$

I am trying to solve the following system of non-linear PDEs with 1 time + 2 space dimensions:

 kp = 0.504; kt = 1.31; kd = 2.8*10^-7;

 equ1 = D[Ini[t, y, z], t] == -kd*Iabs[y, z]*Ini[t, y, z];
 equ2 = D[Rad[t, y, z], t] == 2*kd*Iabs[y, z]*Ini[t, y, z] - 2*kt*(Rad[t, y, z])^2;
 equ3 = D[DB[t, y, z], t] == -kp*Rad[t, y, z]*DB[t, y, z];
 sys = {equ1, equ2, equ3}; 
 iniconditions = {DB[0, y, z] == 12.89*10^3, Rad[0, y, z] == 0, Ini[0, y, z] == 66.2};

The “perturbation” Iabs[x,y] of the system is symmetric:

Input perturbation of the system

Here is the code to generate it:

wo = 5.1*10^-7; zR = 2.49*10^-6; I1 = 9.4*10^6; Eps = 387.6; zo = 200*10^-6;
w[z_] = wo*Sqrt[1 + (z/zR)^2];
Iabs[y_, z_] = I1*(wo/w[z])*Exp[(-2*y^2)/w[z]^2]*Exp[-Eps*(-z + zo)]; 

And I am therefore expecting a symmetric result. NDSolve is able to provide a numeric solution to the system of equation using the "TimeIntegration" Method. However the resulting interpolations are asymmetric:

 sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, {y, -100*wo, 100*wo}, {z, -100*10^-6, 200*10^-6}]
 DB[t_, y_, z_] = DB[t, y, z] /. sol;
 DensityPlot[(1 - DB[0.1, y, z]/(12.89*10^3))*100, {y, -100*wo,100*wo}, {z, -100*10^-6, 100*10^-6}, AspectRatio -> Automatic, PlotLegends -> Automatic, PlotRange -> All,ColorFunction -> Hue, PlotPoints -> 100];

Result of NDSolve

Following this reference I tried to increase the number of plotted points with PlotPoints, it only smoothed out the result. Therefore I am also trying to refine the spatial mesh over the yz domain.

I have not found an option to do it with the TimeIntegration Method and using a rectangle domain to solve the system:

Omega = Rectangle[{-100*wo, -100*10^-6}, {100*wo, 200*10^-6}];
sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, Element[{y,z},Omega], Method -> {"TimeIntegration" -> {"Adams"}}]

Yields the following error messages

NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve

My question is then:

  • Is there a way to controllably semi-discretize the problem with a more refined mesh of the yz spatial domain while maintaining the TimeIntegration method?

Edit: Using ToElementMesh I manually defined the spatial mesh over the yz-domain and passed it to NDSolve:

Needs["NDSolve`FEM`"]
Omega = ImplicitRegion[True, {{y, -100*wo, 100*wo}, {z,-100*10^-6, 200*10^-6}}];
mesh = ToElementMesh[Omega, MaxCellMeasure -> wo/10 000];
mesh["Wireframe"]

With MaxCellMeasure -> wo/10 000, the mesh is a bit coarse but I also tried with MaxCellMeasure-> wo/100 000 unsuccessfully:

When I pass this mesh to NDSolve, I obtain the same error message as when using ImplicitRegion

sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, Element[{y, z}, mesh], Method -> {"TimeIntegration" -> {"Adams"}}];
NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve. >>

So I don't understand why NDSolve is able to solve the system by itself using the TimeIntegration Method but whenever I try to improve the mesh quality it doesn't work.

$\endgroup$
  • $\begingroup$ Something is going wrong with the brackets, did you mean I1*(wo/w[z])*Exp[(-2*y^2)/w[z]^2] Exp[-Eps*(-z + zo)]? $\endgroup$ – Feyre Aug 15 '16 at 8:10
  • $\begingroup$ Yes, thanks for the correction. $\endgroup$ – PAppMath Aug 15 '16 at 8:16
  • $\begingroup$ What code do you actually use to generate the (asymmetric) plot? The NDSolve[] generates three 3D InterpolationFunctions $\endgroup$ – Feyre Aug 15 '16 at 8:48
  • $\begingroup$ Please chek the edited code above the picture of the solution $\endgroup$ – PAppMath Aug 15 '16 at 8:56
  • $\begingroup$ Have you seen ToElementMesh and the related FEM tutorials? $\endgroup$ – Michael E2 Aug 15 '16 at 11:15
3
$\begingroup$

Thanks to @MichaelE2's comment, I tried the following:

sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, {y,-100*wo, 100*wo}, {z, -100*10^-6, 200*10^-6},Method -> {"PDEDiscretization" ->{"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 400}}}]

Which gives a correct symmetric solution to the system of PDEs:

DB[t_, y_, z_] = DB[t, y, z] /. sol;
DensityPlot[(1 - DB[0.3, y, z]/(12.89*10^3))*100, {y, -40*wo,40*wo}, {z,-74*10^-6, 74*10^-6}, AspectRatio -> Automatic, PlotLegends -> Automatic, PlotRange -> All, ColorFunction -> Hue, PlotPoints -> 100]

Correct symmetric solution

The trick was here to put "MinPoints" to at least 200 otherwise NDSolve was returning:

Power::infy: Infinite expression 1/0. encountered. >>
Encountered non-numerical value for a derivative at t == 0.`. >> 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.