0
$\begingroup$

I'm experimenting with different schemes for numerically solving PDEs in Mathematica. I use NestList to implement the Euler Method for 1-D PDE, the code is concise:

forwardEuler[t0_, y0_, tend_, h_, func_] := 
Module[{steps = Floor[(tend - t0)/h], ics = {t0, y0}},
NestList[{#[[1]] + h, #[[2]] + h *func /. {t -> #[[1]], y -> #[[2]]}} &,
ics, steps]]

and it works well.

ListLinePlot[{Table[Evaluate[{t,DSolveValue[{y'[t] == (y[t] + t^2 - 2)/(t + 1), y[0] == 2}, y[t],t]}], {t, 0, 6, 0.2}], 
forwardEuler[0, 2, 6, 0.2, (y + t^2 - 2)/(t + 1)]},PlotStyle -> {Blue, {Dashed, Orange}}, PlotLegends -> {"Exact Solution", "Forward Euler"}]

enter image description here

However, now I want to use Euler Method and central finite difference to solve PDE like the following.

$$ \left\{\begin{array}{l} \frac{\partial u}{\partial t}-\frac{\partial^{2} u}{\partial x^{2}}=x \mathrm{e}^{t}-6 x, \quad 0<x<1, \quad 0<t \leqslant 1, \\ u(x, 0)=x^{3}+x, \quad 0 \leqslant x \leqslant 1, \\ u(0, t)=0, \quad u(1, t)=1+\mathrm{e}^{t}, \quad 0<t \leqslant 1 . \end{array}\right. $$

Is it still recommended to write an Euler method using NestList (or any function from the Nest family)? It has a simple form in 1-D, but I'd like to see it in 2-D or even 3-D.

Additional: I also write codes of Implicit Euler, Trapezoidal, and Heun for 1D equations:

backwardEuler[t0_, y0_, tend_, h_, func_, tol_] := Module[{func1},
  func1[{\[FormalT]_, \[FormalY]_}] := {\[FormalT] + h, NestWhile[N[\[FormalY] + h*func /. {t -> \[FormalT] + h, y -> #}] &, \[FormalY] + h*func /. {t -> \[FormalT], y -> \[FormalY]}, (Abs[#1 - #2] > tol) &, 2, 10]};
  NestList[func1[#] &, {t0, y0}, Floor[(tend - t0)/h]]]

trapezoidal[t0_,y0_,tend_,h_,func_,tol_]:=Module[{func1},
func1[{\[FormalT]_,\[FormalY]_}]:={\[FormalT]+h,NestWhile[N[\[FormalY]+h*((func/.{t->\[FormalT],y->\[FormalY]})+(func/.{t->\[FormalT]+h,y->#}))/2]&,\[FormalY]+h*func/.{t->\[FormalT],y->\[FormalY]},(Abs[#1-#2]>tol)&,2,10]};
NestList[func1[#]&,{t0,y0},Floor[(tend-t0)/h]]]

HuenEuler[t0_,y0_,tend_,h_,func_]:=Module[{func1},
func1[{\[FormalT]_,\[FormalY]_}]:={\[FormalT]+h,N[\[FormalY]+h*((func/.{t->\[FormalT],y->\[FormalY]})+(func/.{t->\[FormalT]+h,y->(\[FormalY]+h*func/.{t->\[FormalT],y->\[FormalY]})}))/2]};
NestList[func1[#]&,{t0,y0},Floor[(tend-t0)/h]]]
Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Euler method etc. are methods that can be extended to work on ODE system naturally, but your implementation can not. What you really need to think about is: how to code Euler method etc. so they work on ODE system? ( A single ODE is just a special case of ODE system. ) Related: mathematica.stackexchange.com/a/274162/1871 $\endgroup$
    – xzczd
    Jan 10 at 10:00

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.