Διαφορικές εξισώσεις

Clear["Global`*"] eq = y'[t] + t^2 y[t] == 1 init = y[0] == 1

Ακριβής επίλυση

sol = DSolve[{eq, init}, y[t], t] // Flatten ySol[t_] := Evaluate[y[t] /. sol] ySol[t] Plot[ySol[t], {t, 0, 10}]

Αριθμητική επίλυση

Nsol = NDSolve[{eq, init}, y[t], {t, 0, 10}]; yNsol[t_] := Evaluate[y[t] /. Nsol] yNsol[2] Plot[yNsol[t], {t, 0, 10}]

Συναρτήσεις Green

Στήσιμο Σ.Δ.Ε. Clear["Global`*"] ODE = -D[p[t] *u'[t], t] + q[t] *u[t] == f[t] init1 = a1 u[a] + a2 u'[a] == 0 init2 = b1 u[b] + b2 u'[b] == 0 (*Πρέπει p>0*) p[t_] := 1 q[t_] := -2 a = 0; b = 1; a1 = 1; a2 = 0; b1 = 1; b2 = 0; ODE init1 init2 GreenFunction[{ODE[[1]], init1, init2}, u[t], {t, a, b}, ξ] // Simplify g[t_, ξ_] := Evaluate[%] g[t, ξ] f[t_] := t^2 ySol[t_] := Integrate[f[ξ]*g[t, ξ], {ξ, a, b}] ySol[t] Απλοποίηση λύσης Refine[ySol[t], t > a && t < b] // Simplify Σύγκριση με ετοιματζίδικη λύση DSolve[{ODE, init1, init2}, u[t], t] // Simplify 1/2 (1 - t^2 + Csc[Sqrt[2]] Sin[Sqrt[2] (-1 + t)]) - 1/2 (1 - t^2 - Cos[Sqrt[2] t] + Cot[Sqrt[2]] Sin[Sqrt[2] t]) // FullSimplify

Ανάπτυγμα σε ιδιοσυναρτήσεις

Ακριβείς ιδιοτιμές/ιδιοσυναρτήσεις

Clear["Global`*"] ODE = a[x]*y''[x] + b[x]*y'[x] + c[x]*y[x] == f[x] Εξειδικεύουμε τις συναρτήσεις a, b, c και το σύνορο L. Ακολούθως προσδιορίζουμε τις πρώτες $n_0$ ιδιοτιμές/ιδιοσυναρτήσεις του διαφορικού τελεστή (T) που συστήνει την Σ.Δ.Ε. μας (Ty=f). a[x_] := 2 b[x_] := -1 c[x_] := 0 L = 2 Pi; n0 = 4; {eigV, eigF} = DEigensystem[{ODE[[1]], DirichletCondition[y[x] == 0, True]}, y[x], {x, 0, L}, n0] Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_n$ οι ιδιοτιμές και $g_n$ οι ιδιοσυναρτήσεις, τότε θέλουμε να βρούμε συντελεστές $y_n$, ώστε: $y(x)=\sum_{n=1}^{\infty} y_n g_n(x)$ Όμως: $Ty(x)=f (x)\Leftrightarrow T(\sum_{n=1}^{\infty} y_n g_n(x))=f (x)\Leftrightarrow \sum_{n=1}^{\infty} y_n Tg_n(x)=f(x)\Leftrightarrow \sum_{n=1}^{\infty} y_n \lambda_n g_n(x)=f(x)$ Σημείο κλειδί στα επόμενα είναι η ορθογωνιότητα των ιδιοσυναρτήσεων, βάσει του εσωτερικού γινομένου: $ =\int_{0}^{L} w(x) f(x) g(x) dx$, όπου $w(x)=\pm \frac{1}{a(x)} \exp(\int \frac{b(x)}{a(x)}dx)$ η (θετική οπωσδήποτε) συνάρτηση βάρους. Για να προσδιοριστούν οι συντελεστές $y_n$ δεν έχουμε παρά να πολλαπλασιάσουμε εσωτερικά την τελευταία ισότητα με $g_n$. Εν προκειμένω θα εξειδικεύσουμε την f. w[x_] := 1/Abs[a[x]] Exp[Integrate[b[x]/a[x], x]] w[x] f[x_] := x^2 yn = Table[Integrate[w[x]*f[x]*Evaluate[eigF[[n]]], {x, 0, L}]/(eigV[[n]]*Integrate[w[x]*Evaluate[eigF[[n]]]^2, {x, 0, L}]), {n, 1, n0}] Η προσεγγιστική μας λύση είναι: yAp[x_] := Evaluate[Sum[yn[[j]]*eigF[[j]], {j, 1, n0}] ] yAp[x] yAp[2] // N Θα την αντιπαραβάλλουμε με τη λύση του `NDSolve`. sol = NDSolve[{ODE, DirichletCondition[y[x] == 0, True]}, y[x], {x, 0, L}]; ySol[x_] := y[x] /. sol[[1]] Plot[{ySol[x], yAp[x]}, {x, 0, L}, PlotLegends -> Automatic]

Προσεγγιστικές ιδιοτιμές/ιδιοσυναρτήσεις

Clear["Global`*"] ODE = a[x]*y''[x] + b[x]*y'[x] + c[x]*y[x] == f[x] init1 = y[0] == 0 init2 = y[L] == 0 Εξειδικεύουμε τις συναρτήσεις a, b, c και το σύνορο L. a[x_] := 2 b[x_] := x c[x_] := 3 L = 2 Pi; Ακολούθως προσδιορίζουμε τις πρώτες $n_0$ ιδιοτιμές/ιδιοσυναρτήσεις του διαφορικού τελεστή (T) που συστήνει την Σ.Δ.Ε. μας (Ty=f). n0 = 30; {eigV, eigF} = NDEigensystem[{ODE[[1]], init1, init2}, y[x], {x, 0, L}, n0]; eigV Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_n$ οι ιδιοτιμές και $g_n$ οι ιδιοσυναρτήσεις, τότε θέλουμε να βρούμε συντελεστές $y_n$, ώστε: $y(x)=\sum_{n=1}^{\infty} y_n g_n(x)$ Όμως: $Ty(x)=f (x)\Leftrightarrow T(\sum_{n=1}^{\infty} y_n g_n(x))=f (x)\Leftrightarrow \sum_{n=1}^{\infty} y_n Tg_n(x)=f(x)\Leftrightarrow \sum_{n=1}^{\infty} y_n \lambda_n g_n(x)=f(x)$ Σημείο κλειδί στα επόμενα είναι η ορθογωνιότητα των ιδιοσυναρτήσεων, βάσει του εσωτερικού γινομένου: $ =\int_{0}^{L} w(x) f(x) g(x) dx$, όπου $w(x)=\pm \frac{1}{a(x)} \exp(\int \frac{b(x)}{a(x)}dx)$ η (θετική οπωσδήποτε) συνάρτηση βάρους. Για να προσδιοριστούν οι συντελεστές $y_n$ δεν έχουμε παρά να πολλαπλασιάσουμε εσωτερικά την τελευταία ισότητα με $g_n$. Εν προκειμένω θα εξειδικεύσουμε την f. w[x_] := Evaluate[1/Abs[a[x]] Exp[Integrate[b[x]/a[x], x]]] w[x] w[2] f[x_] := 2 yn = Table[NIntegrate[w[x]*f[x]*Evaluate[eigF[[n]]], {x, 0, L}]/(eigV[[n]]*NIntegrate[w[x]*Evaluate[eigF[[n]]]^2, {x, 0, L}]), {n, 1, n0}] Η προσεγγιστική μας λύση είναι η: yAp[x_] := Evaluate[Sum[yn[[j]]*eigF[[j]], {j, 1, n0}] ] yAp[2] yn[2] Θα την αντιπαραβάλλουμε με τη λύση του `NDSolve`. sol = NDSolve[{ODE, init1, init2}, y[x], {x, 0, L}]; ySol[x_] := y[x] /. sol[[1]] Plot[yAp[x], {x, 0, L}, PlotLegends -> Automatic] Plot[ySol[x], {x, 0, L}] Plot[{ySol[x], yAp[x]}, {x, 0, L}, PlotLegends -> "Expressions"]

Χωρίς συνάρυηση βάρους

Clear["Global`*"] ODE1 = a[x]*y''[x] + b[x]*y'[x] + c[x]*y[x] == f[x] init1 = y[0] == 0 init2 = y[L] == 0 Μετασχηματίζουμε την Σ.Δ.Ε., με τρόπο που να γίνει Sturm-Liouville. a[x_] := 2 b[x_] := x c[x_] := 3 L = 2 Pi; (*SLoperator = ODE1[[1]]/a[x] Exp[Integrate[b[x]/a[x],x]]*) newF[x_] := Exp[Integrate[b[t]/a[t], {t, 0, x}]] f[x]/a[x] ODE = Exp[Integrate[b[x]/a[x], x]] y''[x] + Exp[Integrate[b[x]/a[x], x]] b[x]/a[x]*y'[x] + Exp[Integrate[b[x]/a[x], x]] c[x]/a[x]*y[x] == Exp[Integrate[b[x]/a[x], x]] f[x]/a[x] n0 = 3; {eigV, eigF} = NDEigensystem[{ODE[[1]], DirichletCondition[y[x] == 0, True]}, y, {x, 0, L}, n0]; (*Simplify[ODE[[1]]/Exp[Integrate[b[x]/a[x],x]]]*) DEigensystem[{ODE[[1]], DirichletCondition[y[x] == 0, True]}, y, {x, 0, L}, n0]; eigF[[2]][5] Plot[eigF[[2]][x], {x, 0, L}] Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_n$ οι ιδιοτιμές και $g_n$ οι ιδιοσυναρτήσεις, τότε θέλουμε να βρούμε συντελεστές $y_n$, ώστε: $y(x)=\sum_{n=1}^{\infty} y_n g_n(x)$. Όμως: $Ty(x)=\frac{f (x)\exp(\frac{b(x)}{a(x)})}{a(x)}\Leftrightarrow T(\sum_{n=1}^{\infty} y_n g_n(x))=\frac{f (x)\exp(\frac{b(x)}{a(x)})}{a(x)}\Leftrightarrow \sum_{n=1}^{\infty} y_n Tg_n(x)=\frac{f (x)\exp(\frac{b(x)}{a(x)})}{a(x)}\Leftrightarrow \sum_{n=1}^{\infty} y_n \lambda_n g_n(x)=\frac{f (x)\exp(\frac{b(x)}{a(x)})}{a(x)}$ Σημείο κλειδί στα επόμενα είναι η ορθογωνιότητα των ιδιοσυναρτήσεων, βάσει του εσωτερικού γινομένου: $ =\int_{0}^{L} f(x) g(x) dx$ Για να προσδιοριστούν οι συντελεστές $y_n$ δεν έχουμε παρά να πολλαπλασιάσουμε εσωτερικά την τελευταία ισότητα με $g_n$. Εν προκειμένω θα εξειδικεύσουμε την f. f[x_] := x^2 (*(f[x]*Exp[Integrate[b[x]/a[x],x]]*eigF[[n]][x])/a[x]*) yn[n_] := NIntegrate[ Exp[Integrate[b[x]/a[x], x]] f[x]/a[x]*eigF[[n]][x], {x, 0, L}]/(eigV[[n]]*NIntegrate[eigF[[n]][x]^2, {x, 0, L}]) For[i = 1, i <= n0, i++, Print[i, ")", yn[i]]] yn[2] Η προσεγγιστική μας λύση είναι: yAp[x_] := Sum[yn[j]*eigF[[j]][x], {j, 1, n0}] yAp[2] yn[2] Θα την αντιπαραβάλουμε με τη λύση του `NDSolve`. sol = NDSolve[{ODE1, init1, init2}, y[x], {x, 0, L}]; ySol[x_] := y[x] /. sol[[1]] Plot[yAp[x], {x, 0, L}, PlotLabel -> "yAp"] Plot[ySol[x], {x, 0, L}, PlotLabel -> "ySol"] Plot[{ySol[x], yAp[x]}, {x, 0, L}, PlotLegends -> "Expressions"]