Método de diferencias finitas para la ecuación de Poisson 1D

Hola usuario de nuestro sitio, hemos encontrado la respuesta a tu pregunta, continúa leyendo y la hallarás aquí.

Solución:

Lo resolví usando FDM y la respuesta es correcta.

ClearAll[u, x];
h = 1/3;
eq1 = -2 u0 + 2 u1 == 0;
eq2 = u0 - 2 u1 + u2 == (6*h)*h^2;
eq3 =  u1 - 2 u2 == (6*2 h - 1/h^2)*h^2;
pts = Solve[eq1, eq2, eq3, u0, u1, u2]

sol = u[x] /. [email protected]





[u''[x] == 6*x, u'[0] == 0, u[1] == 1, u[x], x]
p1 = Plot[sol, x, 0, 1];
p2 = ListPlot[0, 1/9, h, 1/9, 2 h, 1/3, 1, 1,  PlotStyle -> Red];
Show[p1, p2]

El error es grande, ya que $ h $ es largo. Con más puntos, mejorará.

---

Aquí hay un truco rápido para mostrar el efecto de agregar más puntos a FDM

makeA[n_] := Module[A, i, j,
   A = Table[0, i, n, j, n];
   Do[
    Do[
     A[[i, j]] = If[i == j, -2, If[i == j + 1 || i == j - 1, 1, 0]],
     j, 1, n
     ],
    i, 1, n
    ];
   A[[1, 2]] = 2;
   A
   ];
makeB[n_, h_, f_] := Module[b, i,
   b = Table[0, i, n];
   Do[
    b[[i]] = If[i == 1, 0,
      If[i < n, f[(i - 1)*h]*h^2, (f[(i - 1)*h] - 1/h^2)*h^2]
      ]
    , i, 1, n
    ];
   b
   ];
f[x_] := 6*x;(*RHS of ode*)
Manipulate[
 Module[h, A, b, sol, solN, p1, p2, x,
  h = 1/(nPoints - 1);
  A = makeA[nPoints - 1];
  b = makeB[nPoints - 1, h, f];
  sol = LinearSolve[A, b];
  solN = Table[n*h, sol[[n + 1]], n, 0, nPoints - 2];
  AppendTo[solN, 1, 1];
  p1 = Plot[x^3, x, 0, 1];
  p2 = ListLinePlot[solN, PlotStyle -> Red, Mesh -> All];
  Grid[
    Row[" h = ", NumberForm[[email protected], 5, 4]],
    Show[p1, p2, 
      PlotLabel -> "Red is numerical, Blue is exact solution",
      GridLines -> Automatic,
      GridLinesStyle -> LightGray, ImageSize -> 400]
    ]
  ],
 nPoints, 3, "How many points?", 3, 20, 1, Appearance -> "Labeled",
 TrackedSymbols :> nPoints
 ]

tengo 1/9, 1/9, 1/3 lo que está mal debería ser 0, 1/27, 8/27, ya que la exacta $ u = x ^ 3 $.

Como ilustra Nasser, 1/9, 1/9, 1/3 no está mal. Si quieres obtener 0, 1/27, 8/27, una posible solución (no estoy seguro si es la única solución) es usar una fórmula de diferencia de orden superior para aproximar el bc en $ x = 0 $. La fórmula que usaste en $ x = 0 $ es la fórmula de diferencia hacia adelante de primer orden:

$$ f ‘(0) approx frac f (h) -f (0) h $$

Si usamos la fórmula de diferencia unilateral de tercer orden en su lugar:

$$ f ‘(0) approx- frac 11 f (0) 6 h + frac 3 f (h) h – frac 3 f (2 h) 2 h + frac f (3 h) 3 h $$

[email protected]-2 + 9 u[0] - 18 u[1/3] + 9 u[2/3] == 0, 
       -4 + 9 u[1/3] - 18 u[2/3] + 9 u[1] == 0, 
       -((11 u[0])/2) + 9 u[1/3] - 9/2 u[2/3] + u[1] == 0, 
       -1 + u[1] == 0
(*
u[0] -> 0, u[1/3] -> 1/27, u[2/3] -> 8/27, u[1] -> 1
 *)

Por cierto, estas ecuaciones en diferencias se pueden generar fácilmente con mi pdetoae:

eq = u''[x] == 6 x;
bc = u' [0] == 0, u [1] == 1;

grid = Range[0, 1, 1/3];
difforder = 3;

(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[u[x], grid, difforder];

ae = ptoafunc[eq][[2 ;; -2]]
aebc = [email protected][email protected]@ae, aebc

Si quieres el $ A $ y $ b $, Solo usa CoefficientArrays:

barray, marray = CoefficientArrays[[email protected]ae, aebc, u /@ grid];
MatrixForm /@ barray, marray

LinearSolve[marray, -barray]
(* 0, 1/27, 8/27, 1 *)

Con la ayuda de pdetoae, tratar con bcs en (b) y (c) es trivial, por lo que me gustaría omitir el código aquí.

Puedes secundar nuestra función añadiendo un comentario o dejando una puntuación te lo agradecemos.

---

---

Utiliza Nuestro Buscador