Közönséges differenciálegyenletek kezelése Mathematicával

Előkészületek

<< Graphics`

<< Graphics`Colors`

Szimbolikus megoldás

egy = y'[x] == 2  x - y[x]; kezd1 = y[3] == -3; 
    kezd2 = y[2] == -3;

DSolve[{egy, kezd1}, y[x], x]

DSolve[{egy, kezd2}, y[x], x]

ClearAll[ds];
ds[c_] := 
      Simplify[DSolve[{egy, y[1] == c}, y[x], 
            x][[1, 1, 2]]]

ds[c]

DSolve[{y'[x] == a x + b y[x]}, y[x], x]

Iránymező, irányvonalak

pvf = 
        PlotVectorField[{1, 2  x - y}, {x, -2, 4}, {y, -2, 4}, 
          ScaleFactor -> 0.5];

l[k_] := 2  x - k;
iranyok = 
        Plot[{l[-3], l[0], l[1], l[2]}, {x, -2, 4}, 
          PlotRange -> {-2, 4}, 
          PlotStyle -> {{Thickness[0.008], Black}, {Thickness[0.008], 
                Red}, {Thickness[0.008], Black}, {Thickness[0.008], 
                Black}}];

pl[c_] := 
      Plot[Evaluate[ds[c]], {x, -2, 4}, 
        PlotStyle -> {Thickness[
              0.012], ((1 + Sign[c])/2) Green + ((1 - Sign[c])/2) 
                Blue}, PlotRange -> {-2, 4}, 
        DisplayFunction -> Identity]

Show[{pl[2], pl[-3]}, 
        DisplayFunction -> $DisplayFunction];

Show[pvf, iranyok, pl[1], pl[2], pl[3], pl[-1], pl[-2], 
        pl[-4], PlotRange -> {-2, 4}, 
        DisplayFunction -> $DisplayFunction, Frame -> True, 
        Axes -> True];

?*Field*

Elsőrendű lineáris egyenletek

Fönt ilyen szerepelt.

DSolve[{y'[x] + (E^((x^2))) y[x] == 
          Sin[Log[x]], y[1] == 3}, y[x], x]

nds = 
      NDSolve[{y'[x] + (E^((x^2))) y[x] == 
            Sin[Log[x]], y[1] == 3}, y[x], {x, 1, 4}]

Plot[y[x] /. nds, {x, 1, 4}, 
        PlotStyle -> {Dashing[{0.1, 0.02}], Magenta, 
            AbsoluteThickness[3]}, PlotRange -> All];

Az állandók variálásának módszere végigcsinálható.

Egzakt egyenletek

Ezt akarjuk megoldani:

ex = (4  (x^3) y[x]^3 - 
                2  x y[x]) + (3  (x^4) y[x]^2 - (x^2)) y'[
                x] == 0;

P[x_, y_] := 
        ex[[1, {1, 2}]] /. {y[x] -> 
              y};
Q[x_, 
        y_] := ((ex[[1, 3]] /. {y[x] -> y})/(y'))[x]

Keresztbevett parciálisok (+tartomány!!)

D[P[x, y], y] == D[Q[x, y], x]

A primitív függvény meghatározása

Amegoldás, F[x,y]=állandó, szintvonalak

F[x_, y_] = 
      Integrate[Q[x, y], y] + 
        Integrate[P[x, y] - Integrate[D[Q[x, y],  x],  y], x]

Az (1, 1) ponton átmenő megoldások:

Solve[F[x, y] == 0, y]

Solve[F[x, y] == 0, x]

F[x, y]

Plot3D[
        Evaluate[F[x, y]],  {x,  -4, 4},  {y, -1, 
          1}];

ContourPlot[
        F[x,  y], {x, -5, 5}, {y, -5, 5}, Contours -> 20,  
        ContourShading ->  False];

ContourPlot[F[x,  y], {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
        Contours -> 20,  ContourShading ->  False, 
        PlotPoints -> 50];

ContourPlot[F[x,  y], {x, -1.5, 0}, {y, 0, 1.5}, 
        Contours -> 20,  ContourShading ->  False, 
        PlotPoints -> 50];

Integráló tényező

Lehetne ezt is.

Közelítő módszerek

A pontos megoldás

pontos = 
      DSolve[{y'[x] == x y[x],  y[0] == 1}, y[x], 
        x]

Euler-módszer (vö. általános tételek)

euler[f_, h_, x0_, y0_, n_:  3] := 
      Module[{lep, x, y}, lep[{x_, y_}] := {x + h, y + h f[x, y]}; 
        NestList[lep, {x0, y0}, n]]

ClearAll[F];
F[x_, y_] := x y;

lp = 
        ListPlot[euler[F, 0.01, 0, 1, 100], PlotJoined -> True, 
          PlotStyle -> {RGBColor[1, 0, 0]}, 
          DisplayFunction -> Identity];

pl = 
        Plot[(E^(((x^2)/2))), {x, 0, 1}, 
          DisplayFunction -> Identity];

Show[pl, lp, DisplayFunction -> $DisplayFunction];

abra[n_] := 
      Module[{lp, pl}, 
        lp = ListPlot[euler[F, (1. /n), 0, 1, n], PlotJoined -> True, 
            PlotStyle -> {RGBColor[1, 0, 0]}, 
            DisplayFunction -> Identity]; 
pl = 
          Plot[(E^(((x^2)/2))), {x, 0, 1}, 
            DisplayFunction -> Identity]; 
Show[pl, 
          lp, DisplayFunction -> $DisplayFunction, 
          PlotLabel -> "Lépésszám: " <> 
              ToString[n]];]

{abra[5], abra[100]};

Taylorsor-módszer

ts = 
      Series[(E^(((x^2)/2))), {x, 0 , 10}] // 
        Normal

Hatványsormódszer

Szukcesszív approximáció (vö. általános tételek)

Az eszköz

NestList[φ, (Subscript[x, 0]), 3]

nl = NestList[Cos, -6.7, 12]

ListPlot[nl, PlotStyle -> {PointSize[0.02]}];

Alkalmazzuk a differenciálegyenletre

Az egyenértékű integrálegyenlet:

(∫_0^x) 
                y'[ξ] \[DifferentialD]ξ == (∫_0^x) ξ y[ξ] \[DifferentialD]ξ /. {y[0] -> 1}

Az operátor, amit ismételgetünk:

A[φ_] := 1 + Integrate[x φ,  x]

nl = NestList[A ,  1, 2]

Összehasonlítás

Plot[
        Evaluate[{y[x] /. pontos, 
            nl[[3]]}], {x, 0, 4}, 
        PlotStyle -> {Red, Beige}];

Elsőrendű lineáris egyenletrendszerek

DSolve[{x'[t] == -x[t] + 8  y[t], y'[t] == 
                x[t] + y[t]}, {x[t], y[t]}, 
            t][[1]] // Simplify // 
      TableForm

MatrixForm[A = {{-1, 8}, {1, 1}}]

Eigenvalues[A]

Eigenvectors[A]

(E^(Eigenvalues[A] t))

Eigenvectors[A] {(Subscript[d, 1]), (Subscript[d, 2])}

(E^(Eigenvalues[A] t)) . (Eigenvectors[A] {(Subscript[d, 1]), 
            (Subscript[d, 2])})

Ugyanaz, csak másképp vannak jelölve az együtthatók.

Magasabb rendű lineáris egyenletek

.Az átviteli elv programozható. * Egy változó együtthatósat mutatok, ez egy Euler-egyenlet.

DSolve[{(x^2) y''[x] - 4  x y'[x] + 6  y[x] == 
          Log[x], y[1] == 1,  y'[1] == 0}, y[x], 
      x]

Ez meg egy peremértékfeladat. ,,Csoda'', hogy a megoldása egyértelmű!

DSolve[{y''[x] - 2  y'[x] + 2  y[x] == 0, 
        y[0] == 0, y'[π] == (E^π)}, 
      y[x], x]

A Laplace-transzformáció

LaplaceTransform[x'[t], t, 
        p] /. {LaplaceTransform[x[t], t, p] -> X[p]}

LaplaceTransform[x''[t], t, 
        p] /. {LaplaceTransform[x[t], t, p] -> X[p]}

LaplaceTransform[(∫_0^t) f[s] g[t - s] \[DifferentialD]s, 
      t, p]

Alkalmazása lineáris egyenletekre

LaplaceTransform[egy, x, 
        p] /. {LaplaceTransform[y[x], x, p] -> Y[p]}

Solve[%, Y[p]] /. {y[0] -> 1}

InverseLaplaceTransform[%, p, x][[1, 1, 
      2]]

Mit is csináljunk, ha a kezdeti feltétel nem a nullában van megadva?

Alkalmazása lineáris rendszerekre

lte = 
      LaplaceTransform[{x'[t] == -x[t] + 8  y[t], y'[
                t] == x[t] + y[t]}, t, 
          s] /. {LaplaceTransform[x[t], t, s] -> X[s], 
          LaplaceTransform[y[t], t, s] -> Y[s]}

sol = Solve[lte, {X[s], Y[s]}]

InverseLaplaceTransform[sol, s, t]

Alkalmazása bonyolultabb feladatokra

ClearAll[lte];
lte = 
      LaplaceTransform[y'[x] + 4  y[x] == UnitStep[x], x, 
          s] /. {LaplaceTransform[y[x], x, s] -> Y[s]}

ClearAll[sol];
sol = Solve[lte, Y[s]][[1, 1, 
        2]]

ap = Apart[sol /. y[0] -> 2]

InverseLaplaceTransform[ap, s, x]

DSolve[{y'[x] + 4  y[x] == UnitStep[x], y[0] == 2}, 
        y[x], x] // Simplify

(* 
      az elõzõ feladat pontos megoldása*) 
dsol = 
      DSolve[{y'[x] == x y[x],  y[0] == 1}, y[x], 
        x]

Kell-e tudni megoldani differenciálegyenleteket papíron?

Neptun.