Soluções numéricas de Equações Diferenciais Ordinárias

Inicializações

Txt[sty__] := Text[StyleForm[#1, sty], ##2] & @@ ## &
(Subscript[ti, z])_ := 
      Sequence @@ {FontFamily -> Times, FontSlant -> Italic, 
          FontSize -> z}
(Subscript[tb, z])_ := 
      Sequence @@ {FontFamily -> Times, FontWeight -> Bold, 
          FontSize -> z}

Método de Euler

TrF = TraditionalForm

Como integrar aproximadamente uma equação diferencial ordinária:

d((x(t))/d)t == f(x, t)

com condição inicial

x((Subscript[t, o])) == (Subscript[x, o])

Se soubéssemos a solução x(t)  na vizinhança de (Subscript[x, o]) poderíamos desenvolvê-la em série de Taylor e guardar apenas o termo linear em t

(Subscript[x, 1]) == x[t] ≈ Series[x[t], {t, 0, 3}] /. 
          x[0] -> (Subscript[x, o]) /. 
        t -> Δt /. x'[0] :> 
        HF[(Subscript[(Dt[x, t]), t = 0])]

({{x((Subscript[t, o]) + δt) ≈ (Subscript[x, o])+(Subscript[(((dx)/(\(dt\)))), t = 0]) δt+O(δ(t^2))}, { }, {(Subscript[x, 1]) =  x((Subscript[t, o]) + δt)}, {(Subscript[t, 1]) =  (Subscript[t, o])+ δt   }})

Como (Subscript[(((dx)/(dt))), t = (Subscript[t, o])])=f((Subscript[x, o]),  (Subscript[t, o])), podemos escrever

({{(Subscript[x, 1]) =  
                (Subscript[x, o]) + f((Subscript[x, o]),  (Subscript[t, o])) δt }, {(Subscript[t, 1])= δt  }})

A seguir assumimos que x((Subscript[t, 1]))  ≈  (Subscript[x, 1]) e usamos (Subscript[(((dx)/(dt))), t = (Subscript[t, 1])])=f((Subscript[x, 1]),  (Subscript[t, 1])) para deduzir que

({{x((Subscript[t, 1]) +  δt) ≈ x((Subscript[t, 1]))+(Subscript[(((dx)/(\(dt\)))), t = (Subscript[t, 1])])δt+O(δ(t^2))}, {⇓}, {(Subscript[x, 2]) =  x((Subscript[t, 1]) + δt)}, {(Subscript[t, 2]) =  (Subscript[t, 1])+ δt   }})

pelo que

(Subscript[x, 2])  =  (Subscript[x, 1]) + f((Subscript[x, 1]), (Subscript[t, 1])) δt

Em geral chegamos à expressão

({{(Subscript[x, k + 1]) =  (Subscript[x, k]) + f((Subscript[x, k]), (Subscript[t, k])) δt}, {(Subscript[t, k + 1]) =  (Subscript[t, k]) +  δt}})

Como exemplo vejamos o caso de :

f[x_, t_] := -5 t (x^2)

Uma solução exacta existe para esta equação separável com condição inicial x[0] = 1:

X[t_]=(x[t]/.DSolve[{x'[t]==f[x[t], t],x[0] == 1}//RH, x[t],t])[[1]]

pl = 
        Plot[X[t], {t, -1, 5}, 
          PlotStyle -> {Thickness[0.005`], RGBColor[0, 0, 1]}, 
          PlotPoints -> 100];

A função clr é definida aqui porque é necessário especificar que um subconjunto de definições atribuídas ao símbolo Subscript devem ser Unset. Todas as definições de funções como (Subscript[x, n])_ são atribuídas a Subscript. Fazer Clear[Subscript] apagaria mais do que seria desejável...

\(\(clr := {Off[Unset::"norep"]; 
ClearAll[t, x]; 
x\_n_ =. ; 

          t\_n_ =. ; 
\(HoldForm[Unset[Subscript[x, #]]] &\) /@ Range[0, m] // 
            ReleaseHold; 
\(HoldForm[Unset[Subscript[t, #]]] &\) /@ 
              Range[0, m] // ReleaseHold;};\)\)

Definition[Subscript]

Nas definições que implementam o método de Euler simples, as atribuições da forma (Subscript[x, n])_ :=  (Subscript[x, n]) = def[n] criam uma memória dos valores já calculados para (Subscript[x, k]) à medida que (Subscript[x, n]) é invocado para valores específicos de n = k . Da próxima vez que (Subscript[x, n])for invocado com n = k, def[k] não é recalculado, e o  (Subscript[x, k])memorizado é devolvido.  A seguir (Subscript[lp, n])_ implementa um ListPlot dos valores (Subscript[x, n]) calculados para n instantes (Subscript[t, n]).

clr;

    
δt  :=  .2;
(Subscript[t, 0])  := 0.2;
(Subscript[t, n])_ := (Subscript[t, n]) = (Subscript[t, n - 1]) + δt;
(Subscript[x, 0])  := (Subscript[x, 0]) = X[(Subscript[t, 0])];
\(\(\(x\_n_ := \(x\_n = 
        x\_\(n - 1\) + f[x\_\(n - 1\), t\_\(n - 1\)] δt\)\)\(
\)
    \)\)
\(\(lp\_n_ := 
        ListPlot[\({t\_#, x\_#} &\) /@ Rg[0, n], 
          DisplayFunction -> Identity, Axes -> True, 
          PlotStyle -> {Thickness[0.005`], 
              RGBColor[1, 1, 1]}];\)\)

m = 10;
eul1= SH[{ pl,  (Subscript[lp, m])},
DisplayFunction -> $DisplayFunction,PlotRange -> {{- .8, (m + 2) δt}, {- .1, 
                  1.15}},Axes -> {True, True},
Epilog->{\({AbsoluteDashing[{5, 5, 
                        5}], \(Line[{{t\_#, 0}, {t\_#, x\_#}}] &\) /@ 
                      Range[0, m]}\),RGBColor[1, 0, 0],
Txt[(Subscript[ti, 14])]/@\({Sequence @@ \((\({\(ToString[t]\)\_#, {t\_#, 
                                  0}, {\(-1\), \(-1\)}} &\) /@ 
                            Range[0, 
                              m])\), 
Sequence @@ \((\({\(ToString[x]\)\_#, {t\_#, x\_# +  .01}, {\(-1\), \(-1\)}} &\) /@ 
                            Range[0, m])\)}\)},Background -> GrayLevel[0]];

O método de Euler melhorado

Notando que o integral da equação diferencial entre instantes (Subscript[t, k]) e (Subscript[t, k + 1]) é a área sob a solução verdadeira x(t)

(∫_((Subscript[t, k]))^((Subscript[t, k + 1])))d((x(
                        t))/d)t dt = x((Subscript[t, k + 1]))-x((Subscript[t, k])) =(∫_((Subscript[t, k]))^((Subscript[t, k + 1])))f(x, t) dt

poderíamos dar uma melhor aproximação desta área fazendo

({{ (∫_((Subscript[t, k]))^((Subscript[t, k + 1])))f(x, t) dt ≈ (1/2)(f((Subscript[x, k + 1]), (Subscript[t, k + 1])) + 
                        f((Subscript[x, k]), (Subscript[t, k])))((Subscript[t, k + 1])-(Subscript[t, k]))}, {⇓}, {(Subscript[x, k + 1])= (Subscript[x, k])+(1/2)(f((Subscript[x, k + 1]), (Subscript[t, k + 1])) + f((Subscript[x, k]), (Subscript[t, k]))) δt}, {(Subscript[t, k + 1]) =  (Subscript[t, k])  +  δt}})

Resta agora substituir (Subscript[x, k + 1]) no lado direito por uma primeira aproximação, que pode ser \(TraditionalForm\`x\&~\_\(k + 1\) = \ 
        x\_k + \(f(x\_k, t\_k)\)\ δt\), donde temos

({{\(x\&~\_\(k + 1\) = \ x\_k + \(f(x\_k, t\_k)\)\ δt\)}, { }, {(Subscript[x, k + 1])= (Subscript[x, k])+(1/2)\((f(x\&~\_\(k + 1\), t\_\(k + 1\)) + 
                        f(x\_k, t\_k))\) δt}, { }, {(Subscript[t, k + 1]) =  (Subscript[t, k])  +  δt}})

A implementação do método melhorado é semelhante à anterior:

clr;
δt :=  .2;
(Subscript[t, 0]) := 0.2;
(Subscript[t, n])_ := (Subscript[t, n]) = (Subscript[t, n - 1]) + δt;
(Subscript[x, 0]) := (Subscript[x, 0]) = X[(Subscript[t, 0])];
\(x\_n_ := \(x\_n = 
        x\_\(n - 1\) + \(1\/2\) \((f[
                  x\_\(n - 1\) + f[x\_\(n - 1\), t\_\(n - 1\)] δt, 
                  t\_n] + 
                f[x\_\(n - 1\), 
                  t\_\(n - 1\)])\) δt\)
\)
\(\(lp\_n_ := 
        ListPlot[\({t\_#, x\_#} &\) /@ Rg[0, n], 
          DisplayFunction -> Identity, 
Axes -> 
            True, PlotStyle -> {Thickness[0.005`], 
              RGBColor[1, 1, 0]}];\)\)

m = 10;
eul2=SH[{pl, (Subscript[fp, m])},
DisplayFunction -> $DisplayFunction,PlotRange -> {{- .8, (m + 2) δt}, {- .1, 
                  1.1}},AspectRatio -> 1,Axes -> {True, True},
Epilog->{\({AbsoluteDashing[{5, 5, 
                        5}], \(Line[{{t\_#, 0}, {t\_#, x\_#}}] &\) /@ 
                      Range[0, m]}\),RGBColor[1, 0, 0],
Txt[(Subscript[ti, 14])]/@\({Sequence @@ \((\({\(ToString[t]\)\_#, {t\_#, 
                                  0}, {\(-1\), \(-1\)}} &\) /@ 
                            Range[0, 
                              m])\), 
Sequence @@ \((\({\(ToString[x]\)\_#, {t\_#, x\_# +  .01}, {\(-1\), \(-1\)}} &\) /@ 
                            Range[0, m])\)}\)},Background -> GrayLevel[0]];

A comparação dos dois métodos é feita sobrepondo os gráficos eul1 e eul2.

Show[{eul2, eul1}];

Exemplo: SetDelayed com memória

Clear[g]

g[x_] := g[x] = κ (x^2)

?? g

g[2]

?? g

g[x_] := g[x] = κ (x^2) - 3

?? g

g[3]

g[2]