Recently I asked a question on solving a system of 3 nonlinear PDE (To solve a system of 3 nonlinear PDE). Now I transitioned to the next step, that is, a realistic PDE which has been my initial aim:
eq1A = -0.00006 p1[t, x, y] - 0.00067 p1[t, x, y]^3 +
0.019 p1[t, x, y]^5 - 0.0026 p1[t, x, y] p2[t, x, y]^2 +
0.016 p1[t, x, y]^3 p2[t, x, y]^2 +
0.008 p1[t, x, y] p2[t, x, y]^4 + (y p3[t, x, y])/(x^2 + y^2) -
0.0026 p1[t, x, y] p3[t, x, y]^2 +
0.016 p1[t, x, y]^3 p3[t, x, y]^2 +
0.008` p1[t, x, y] p3[t, x, y]^4 + D[p1[t, x, y], x, x] +
D[p1[t, x, y], y, y] == D[p1[t, x, y], t];
eq2A = -0.00006 p2[t, x, y] - 0.0026 p1[t, x, y]^2 p2[t, x, y] +
0.008 p1[t, x, y]^4 p2[t, x, y] - 0.00068 p2[t, x, y]^3 +
0.016 p1[t, x, y]^2 p2[t, x, y]^3 + 0.02` p2[t, x, y]^5 - (
x p3[t, x, y])/(x^2 + y^2) - 0.0026 p2[t, x, y] p3[t, x, y]^2 -
0.2 p1[t, x, y]^2 p2[t, x, y] p3[t, x, y]^2 +
0.016 p2[t, x, y]^3 p3[t, x, y]^2 +
0.0078 p2[t, x, y] p3[t, x, y]^4 + D[p2[t, x, y], x, x] +
D[p2[t, x, y], y, y] == D[p2[t, x, y], t];
eq3A = (y p1[t, x, y])/(x^2 + y^2) - (x p2[t, x, y])/(x^2 + y^2) -
0.00006 p3[t, x, y] - 0.0026 p1[t, x, y]^2 p3[t, x, y] +
0.008 p1[t, x, y]^4 p3[t, x, y] -
0.0026 p2[t, x, y]^2 p3[t, x, y] -
0.2 p1[t, x, y]^2 p2[t, x, y]^2 p3[t, x, y] +
0.008 p2[t, x, y]^4 p3[t, x, y] - 0.00068 p3[t, x, y]^3 +
0.016 p1[t, x, y]^2 p3[t, x, y]^3 +
0.016 p2[t, x, y]^2 p3[t, x, y]^3 + 0.02 p3[t, x, y]^5 +
D[p3[t, x, y], x, x] + D[p3[t, x, y], y, y] == D[p3[t, x, y], t];
within the disk with the mesh:
Needs["NDSolve`FEM`"];
R = 40.;
bM = ToBoundaryMesh[Disk[{0, 0}, R]];
hCenter = 0.5; (* ~0.1\[Ellipsis]1 near r=0*)
hBoundary = 30; (* ~5\[Ellipsis]10 near r=R*)
h[r_] := hCenter + (hBoundary - hCenter) (r/R)^1.5;
Amax[r_] := Sqrt[3]/4 h[r]^2;
c=Sqrt[3]/4;
mrf = Function[{verts,
area},(*refine if the current area is above the radius-
dependent target*)area > Amax[Norm[Mean[verts]]]];
mesh = ToElementMesh[bM, MeshRefinementFunction -> mrf,
MaxCellMeasure -> c hBoundary^2, "MeshOrder" -> 1];
mesh["Wireframe"]
and boundary as well as initial conditions:
dC1 = DirichletCondition[p1[t, x, y] == 0, ElementMarker == 1];
dC2 = DirichletCondition[p2[t, x, y] == 0, ElementMarker == 1];
dC3 = DirichletCondition[p3[t, x, y] == 0, ElementMarker == 1];
iC1A = p1[0, x, y] == 0.01*Exp[-(x^2 + y^2)/4];
iC2A = p2[0, x, y] == 0.02*Exp[-(x^2 + y^2)/4];
iC3A = p3[0, x, y] == 0.03 Exp[-(x^2 + y^2)/4];
All this worked for me in the previous case. The present case differs from the previous one by the considerably more complex polynomials representing the nonlinear parts of the equations.
Then I apply the standard NDSolve:
ndsl = NDSolve[{eq1A, eq2A, eq3A, dC1, dC2, dC3, iC1A, iC2A,
iC3A}, {p1, p2, p3}, {t, 0, 10}, {x, y} \[Element] mesh];
I get the following response: LinearSolve::mcovl: The computation encountered machine-number overflow. Only machine-number code is available for sparse matrices. If your matrix is not too large, consider trying again after using Normal on the matrix.
Any idea?
D[p1[t, x, y] t],D[p2[t, x, y] t],D[p3[t, x, y] t]. It should be with,asD[p1[t, x, y], t]... $\endgroup$nr = ToNumericalRegion[Disk[{0, 0}, R]]; bM = ToBoundaryMesh[nr]; mesh = ToElementMesh[nr, MeshRefinementFunction -> mrf, MaxCellMeasure -> c hBoundary^2, "MeshOrder" -> 1];This allows the full mesh generation to get access to the symbolic region. I can then generate a better boundary approximation. $\endgroup$cis not defined and why have you chosen to restrict the mesh order to 1? $\endgroup$