8
$\begingroup$

The problem modeled by the Poisson PDE is related to the torsion of prismatic beams and I use the Finite Differences Method (FDM). I've managed to solve the equation over a rectangle region with quadratic mesh, here's my code (the dimensions are 0.2 m × 0.2 m):

a = 0.2; b = 0.2; n = 4; m = 4; deltax = a/n; deltay = b/m; G = 10^6; theta = 0.01745; phi[0, j_] = 0; phi[n, j_] = 0; phi[i_, 0] = 0; phi[i_, m] = 0; fi[i_, j_] = -2 G theta; vars = Flatten[Table[phi[i, j], {i, 1, n - 1}, {j, 1, m - 1}]]; eqns = Flatten[ Table[(phi[i + 1, j] - 2 phi[i, j] + phi[i - 1, j])/ deltax^2 + (phi[i, j + 1] - 2 phi[i, j] + phi[i, j - 1])/ deltay^2 == fi[i, j], {i, 1, n - 1}, {j, 1, m - 1}]]; sol = Solve[eqns, vars][[1]]

But I also need to solve for a circular region with radius 0.1m, and I frankly don't have an idea what to do:

enter image description here

The corresponding difference equation is:

$\frac{2}{\Delta x^{2}}\left(\frac{1}{a(1+a)}\phi_{i-1,j}+\frac{1}{1+a}\phi_{i+1,j}-\frac{1}{a}\phi_{i,j}\right)+\frac{2}{\Delta y^{2}}\left(\frac{1}{b(1+b)}\phi_{i,j-1}+\frac{1}{1+b}\phi_{i,j+1}-\frac{1}{b}\phi_{i,j}\right)=-2G\theta$

I know how to solve it with NDSolve:

R = 0.1; G = 1000000; θ = 0.01745329252; Ω = Disk[{0, 0}, {R, R}]; op = Laplacian[u[x, y], {x, y}] +2 G θ; Subscript[Γ,D] = {DirichletCondition[u[x, y] == 0, True]}; Φ = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y} ∈ Ω]; Plot3D[Φ[x, y], {x, y} ∈ Ω]

enter image description here

I just want to solve the problem with FDM using Cartesian coordinates.

Improve this question
$\endgroup$
8
  • $\begingroup$ No worries :) . $\endgroup$ Commented Oct 31, 2014 at 19:09
  • $\begingroup$ Well, is it necessary for you to use FDM rather than NDSolve solving this equation? To say at least, is it necessary for you to use Cartesian form of the equation rather than polar form? The handling for irregular mesh in FDM is really cumbersome in my view, in fact that's where I stoped my self-learning for FDM. (You may be interested in this post. ) $\endgroup$ Commented Nov 1, 2014 at 3:57
  • $\begingroup$ Yes, because I want to solve the problem with FDM, the solution with NDSolve already I have. $\endgroup$ Commented Nov 1, 2014 at 6:52
  • $\begingroup$ @wlkyr So the polar form isn't what you need either and you simply want to explore FDM for irregular region? $\endgroup$ Commented Nov 1, 2014 at 17:47
  • $\begingroup$ Yes, because I want to solve with Cartesian form, and I find the solution numericali but I dont have idea how to solve with Mathematica $\endgroup$ Commented Nov 1, 2014 at 17:49

1 Answer 1

Reset to default
12
$\begingroup$

As mentioned in the comment above, handling irregular region with FDM is cumbersome and frustrating in my view, actually that's where I stopped my self-learning of FDM and turned to finite element method (FEM), which is more suitable for this task, and finally write this. (Have a look at this post for more information. ) Nevertheless, there seems to be no implementation of FDM on irregular region with Mathematica hitherto, so I'd like to have a try.

Well, I believe that as a centuried technique, there must be some legacy code (probably in other programming language) implementing FDM on irregular region, but strangely I never found a piece. (My skill of googling is too bad?) So the following code is completely self-made and quite Mathematica-style (I think).

The method I chose to treat the irregular boundary is the one OP mentioned in his question: extrapolation. Yeah, that graph and the difference equation illustrates extrapolation. If you're unfamiliar with it, a relatively good learning material is Section 3.4 from p71 and Section 6.4 from p199 of Numerical Solution of Partial Differential Equations - An Introduction (Morton K., Mayers D)… OK, let me explain it briefly.

Graph 1

enter image description here

Red points for the unknowns, gray points for those outside the domain, values on green curve is known.

FDM is a technique solving unknown grids with known grids, but as shown in Graph 1, when implementing FDM on an irregular region, grid points are usually not exactly on the boundary where the function value is known, so what to do? The simplest solution is to treat the grids near the boundary as grids on the boundary, which is called zero-order interpolation and discouraged. (Those material I found doesn't mention the reason, I guess it's probably because the error caused by zero-order interpolation is large and not quantitatively controllable.) To achieve the quantitative control on the error caused by FDM, extrapolation is needed. For example, as illustrated in Graph 2 below, we can use the values on B, C, D to fit a parabola to estimate the value on A:

Graph 2

enter image description here

An equivalent explanation for the extrapolation is weighted difference formula, which is also covered in the book linked above. In the rest part of this answer I'll concentrate on programming. Here's the code:

G = 1000000; θ = Rationalize[0.01745329252, 0]; m = n = 24; R = 1/10; dx = 2 R/m; dy = 2 R/n; r1 = Solve[(x + R)/dx == i, x][[1]]; r2 = Solve[(y + R)/dy == j, y][[1]]; ruleInner = (x^2 + y^2 < R^2 /. r1 /. r2); ruleOuter = (x^2 + y^2 > R^2 /. r1 /. r2); areaInner = Table[ruleInner, {i, 0, m}, {j, 0, n}]; areaOuter = Table[ruleOuter, {i, 0, m}, {j, 0, n}]; ArrayPlot /@ Boole[{areaInner, areaOuter}]

enter image description here

I rationalized all the parameters for convenience. areaInner and areaOuter are masks for extracting the desired grid points:

{var, varNot} = Flatten@Pick[Table[f[i, j], {i, 0, m}, {j, 0, n}], #] & /@ {areaInner, areaOuter};

var are grid points inside the circle i.e. the red points in Graph 1, varNot are the grays.

df[dx_, pos_] := (MapAt[# - 1 &, #, pos] - 2 # + MapAt[# + 1 &, #, pos])/dx^2 & {d2fx, d2fy} = {df[dx, 1] /@ var, df[dy, 2] /@ var};

d2fx and d2fy are $\frac{\partial ^2\phi }{\partial x^2}$ and $\frac{\partial ^2\phi }{\partial y^2}$.

{d2fxOuter, d2fyOuter} = Intersection[varNot, Cases[#, f[__], Infinity]] & /@ {d2fx, d2fy}; {d2fxBorder, d2fyBorder} = Thread[Complement[Cases[#, f[__], Infinity], var, varNot] -> 0] & /@ {d2fx, d2fy};

d2fxOuter and d2fyOuter are grids outside the region but appearing in the equations i.e. As in Graph 2, d2fxBorder and d2fyBorder are those exactly on the boundaries, where As and Bs superpose. (It is relatively rare, actually there're only 4 in this case). Then we use extrapolation to approximate those outside the region:

{pointBorderd2fx, pointBorderd2fy} = With[{ex = #[[1]], pos = #[[2]]}, With[{pos2 = 2 - pos + 1}, {#, 0} & /@ Flatten[Nearest[{i, j}[[pos]] /. Solve[x^2 + y^2 == R^2 /. r1 /. r2, {i, j}[[pos]]] /. {i, j}[[pos2]] -> #[[pos2]], #[[pos]]] & /@ ex]]] & /@ {{d2fxOuter, 1}, {d2fyOuter, 2}};

pointBorderd2fx and pointBorderd2fy are points on the boundary that will be used in extrapolation i.e. Bs in Graph 2.

{pointInnerd2fx, pointInnerd2fy} = With[{ex = #[[1]], pos = #[[2]]}, Map[{#[[pos]], #} &, Intersection[#, var] & /@ Thread /@ MapAt[{-2, -1, 1, 2} + # &, ex, {All, pos}], {2}]] & /@ {{d2fxOuter, 1}, {d2fyOuter, 2}};

pointInnerd2fx and pointInnerd2fy are points inside the domain that will be used in extrapolation i.e. Cs and Ds in Graph 2.

extra[{x4_, y4_}, {x1_, y1_}, {{x2_, y2_}, {x3_, y3_}}] = y4 -> InterpolatingPolynomial[{{x1, y1}, {x2, y2}, {x3, y3}}, x4];

extra provides the extrapolating formula. The last step is to form the equation set with difference equations, eliminating As with Bs, Cs and Ds and then solve it:

sol = First@ NSolve[(d2fx /. d2fxBorder /. MapThread[extra, {{#[[1]], #} & /@ d2fxOuter, pointBorderd2fx, pointInnerd2fx}]) + (d2fy /. d2fyBorder /. MapThread[extra, {{#[[2]], #} & /@ d2fyOuter, pointBorderd2fy, pointInnerd2fy}]) == -2 G θ // Thread, var]; // AbsoluteTiming

Finally let's draw a picture for the solution. Notice that we need to add points on the boundary because they are almost not included in sol:

bound = Flatten /@ Transpose@{#, ConstantArray[0., Length@#]} &@ Cases[Normal@ ContourPlot[ x^2 + y^2 == R^2 /. r1 /. r2 // Evaluate, {i, 0, m}, {j, 0, n}], Line[a__] -> a, ∞][[1]]; ListPlot3D[Join[bound, sol /. (f[a_, b_] -> c_) :> {a, b, c}]]

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ Useful post. Any reference/idea on how to apply extrapolation at a corner node in a 9-point compact stencil discretizing the anisotropic diffusion operator on non-rectangular domain with non-uniform grid? $\endgroup$ Commented Dec 8, 2014 at 18:45
  • $\begingroup$ @unlikely Glad it help. As mentioned, this post is my first trial for FDM on irregular region, and in fact all the references in my hand has been linked in this post, maybe you can consider asking a separate question for your specific equation? $\endgroup$ Commented Dec 11, 2014 at 6:39
  • $\begingroup$ Thanks for help, but I have another question how to solve same example in polar coordinates?please help me $\endgroup$ Commented Mar 28, 2015 at 17:35
  • $\begingroup$ Thanks Jens, I try that in polar coordinates but I have a problem with boundary values. I posted that mathematica.stackexchange.com/questions/78468/… $\endgroup$ Commented Mar 28, 2015 at 21:39
  • $\begingroup$ @xzczd Thanks for the answer, but I have another question with your programs, when I changed m and n the radius of circle changed $\endgroup$ Commented Apr 2, 2015 at 10:58

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.