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

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}]
\[{y(t)\to -\frac{(E^{-\frac{1}{3} (t^{3})}) (-3 (t^{2})+({(-1)}^{2/3}) (3^{1/3}) (t^{2}) Gamma[\frac{1}{3}]-(3^{1/3}) ({(-(t^{3}))}^{2/3}) Gamma[\frac{1}{3},-\frac{1}{3} (t^{3})])}{3 (t^{2})}}\]
\[-\frac{(E^{-\frac{1}{3} (t^{3})}) (-3 (t^{2})+({(-1)}^{2/3}) (3^{1/3}) (t^{2}) Gamma[\frac{1}{3}]-(3^{1/3}) ({(-(t^{3}))}^{2/3}) Gamma[\frac{1}{3},-\frac{1}{3} (t^{3})])}{3 (t^{2})}\]
0 2.5 5 7.5 10 0 0.5 1 1.5

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

Nsol = NDSolve[{eq, init}, y[t], {t, 0, 10}];
yNsol[t_] := Evaluate[y[t] /. Nsol]
yNsol[2]
Plot[yNsol[t], {t, 0, 10}]
{0.4532145795303309`}
0 2.5 5 7.5 10 0 0.5 1 1.5

Συναρτήσεις 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
\[q(t) u(t)-(p')[t] (u')[t]-p(t) (u'')[t]=f(t)\]
\[a1 u(a)+a2 (u')[a]=0\]
\[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
\[-2 u(t)-(u'')[t]=f(t)\]
\(u(0)=0\)
\(u(1)=0\)
GreenFunction[{ODE[[1]], init1, init2}, u[t], {t, a, b},  ξ] // Simplify
\[-\frac{\csc(\sqrt{2}) (HeavisideTheta[-t+\xi ] \sin((\sqrt{2}) t) \sin((\sqrt{2}) (-1+\xi ))+HeavisideTheta[t-\xi ] \sin((\sqrt{2}) (-1+t)) \sin((\sqrt{2}) \xi ))}{\sqrt{2}}\]
g[t_, ξ_] := Evaluate[%]
g[t, ξ]
\[-\frac{\csc(\sqrt{2}) (HeavisideTheta[-t+\xi ] \sin((\sqrt{2}) t) \sin((\sqrt{2}) (-1+\xi ))+HeavisideTheta[t-\xi ] \sin((\sqrt{2}) (-1+t)) \sin((\sqrt{2}) \xi ))}{\sqrt{2}}\]
f[t_] := t^2
ySol[t_] := Integrate[f[ξ]*g[t, ξ], {ξ, a, b}]
ySol[t]
\[\frac{\csc(\sqrt{2}) (-HeavisideTheta[1-t] (-(\sqrt{2}) \cos(\sqrt{2}) HeavisideTheta[-t]-(-1+HeavisideTheta[-t]) ((\sqrt{2}) (-1+t^{2}) \cos((\sqrt{2}) (-1+t))-2 t \sin((\sqrt{2}) (-1+t)))) \sin((\sqrt{2}) t)-HeavisideTheta[t] \sin((\sqrt{2}) (-1+t)) (-(\sqrt{2})+(\sqrt{2}) (-1+t^{2}) \cos((\sqrt{2}) t) (-1+HeavisideTheta[-1+t])+2 HeavisideTheta[-1+t] \sin(\sqrt{2})-2 t (-1+HeavisideTheta[-1+t]) \sin((\sqrt{2}) t)))}{2 (\sqrt{2})}\]

Απλοποίηση λύσης

Refine[ySol[t], t > a && t < b] // Simplify
\[(\frac{1}{2}) (1-t^{2}+\csc(\sqrt{2}) \sin((\sqrt{2}) (-1+t)))\]

Σύγκριση με ετοιματζίδικη λύση

DSolve[{ODE, init1, init2}, u[t], t] // Simplify
\[{{u(t)\to (\frac{1}{2}) (1-t^{2}-\cos((\sqrt{2}) t)+\cot(\sqrt{2}) \sin((\sqrt{2}) t))}}\]
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
0

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

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

Clear["Global`*"]
ODE = a[x]*y''[x] + b[x]*y'[x] + c[x]*y[x] == f[x]
\[c(x) y(x)+b(x) (y')[x]+a(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]
\[{{-\frac{5}{8},-\frac{17}{8},-\frac{37}{8},-\frac{65}{8}},{(E^{x/4}) \sin(\frac{x}{2}),(E^{x/4}) \sin(x),(E^{x/4}) \sin(\frac{3 x}{2}),(E^{x/4}) \cos(x) \sin(x)}}\]

Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_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}]
\[(\frac{1}{2}) (E^{-x/2})\]
\[{-\frac{64 (-32+4 (E^{-\pi /2}) (-8+40 \pi +25 ({\pi }^{2})))}{625 \pi },-\frac{128 (-416+(E^{-\pi /2}) (416-68 \pi (8+17 \pi )))}{83521 \pi },-\frac{192 (-1056+(E^{-\pi /2}) (-1056+148 \pi (8+37 \pi )))}{1874161 \pi },\frac{2048 (E^{-\pi /2}) (-488+488 (E^{\pi /2})+520 \pi +4225 ({\pi }^{2}))}{17850625 \pi }}\]

Η προσεγγιστική μας λύση είναι:

yAp[x_] := Evaluate[Sum[yn[[j]]*eigF[[j]], {j, 1, n0}] ]
yAp[x]
yAp[2] // N
\[-\frac{64 (E^{x/4}) (-32+4 (E^{-\pi /2}) (-8+40 \pi +25 ({\pi }^{2}))) \sin(\frac{x}{2})}{625 \pi }-\frac{128 (E^{x/4}) (-416+(E^{-\pi /2}) (416-68 \pi (8+17 \pi ))) \sin(x)}{83521 \pi }+\frac{2048 (E^{-\frac{\pi }{2}+\frac{x}{4}}) (-488+488 (E^{\pi /2})+520 \pi +4225 ({\pi }^{2})) \cos(x) \sin(x)}{17850625 \pi }-\frac{192 (E^{x/4}) (-1056+(E^{-\pi /2}) (-1056+148 \pi (8+37 \pi ))) \sin(\frac{3 x}{2})}{1874161 \pi }\]
-10.315277653368113`

Θα την αντιπαραβάλλουμε με τη λύση του 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]
0 2 4 6 -20 -10 0
Series 1
Series 2

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

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
\[c(x) y(x)+b(x) (y')[x]+a(x) (y'')[x]=f(x)\]
\(y(0)=0\)
\(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
{0.9970034165921076`,-1.1005977918205117`,-3.6976425762279757`,-7.18772649626623`,-11.673883893355718`,-17.16681356948949`,-23.670949702146633`,-31.19528674854776`,-39.755099154226244`,-49.373039070585705`,-60.08016417074049`,-71.91674151629005`,-84.93248802490166`,-99.18551336010583`,-114.73810693394779`,-131.64395210321698`,-149.90795683103474`,-169.33488564395674`,-188.72726192561353`,-200.14236728467554`,-254.81296842277328`,-281.40179370131506`,-312.8447298044985`,-348.11827243664374`,-387.21456929697706`,-430.35282464418293`,-477.8077043034215`,-529.8348109938067`,-586.6066011056157`,-648.1353528955397`}

Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_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]
\[(\frac{1}{2}) (E^{(\frac{1}{4}) (x^{2})})\]
\[\frac{E}{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}]
{10.972535165512465`,-18.122106200353212`,-7.352530718810228`,-3.8009947199391982`,-2.2100355232407214`,-1.3393665302299378`,-0.8907873585398995`,-0.6053641501120237`,-0.4418196073271723`,-0.3249914161278218`,0.25316242628614744`,0.197340030716863`,-0.16123561207918144`,0.1313889350234795`,-0.11127681497360066`,-0.093985471293157`,0.08208842329194573`,0.07195724067807896`,-0.06803714719077926`,0.04941261369740445`,-0.023666613445375486`,0.02445110286244254`,-0.021710301100581378`,-0.018436566765873196`,-0.01538473767904688`,-0.012590777825525088`,-0.010225952527265082`,0.008212512360652976`,0.006562014679639082`,-0.005211121151164438`}

Η προσεγγιστική μας λύση είναι η:

yAp[x_] := Evaluate[Sum[yn[[j]]*eigF[[j]], {j, 1, n0}] ]
yAp[2]
yn[2]
13.186581524835152`
{10.972535165512465`,-18.122106200353212`,-7.352530718810228`,-3.8009947199391982`,-2.2100355232407214`,-1.3393665302299378`,-0.8907873585398995`,-0.6053641501120237`,-0.4418196073271723`,-0.3249914161278218`,0.25316242628614744`,0.197340030716863`,-0.16123561207918144`,0.1313889350234795`,-0.11127681497360066`,-0.093985471293157`,0.08208842329194573`,0.07195724067807896`,-0.06803714719077926`,0.04941261369740445`,-0.023666613445375486`,0.02445110286244254`,-0.021710301100581378`,-0.018436566765873196`,-0.01538473767904688`,-0.012590777825525088`,-0.010225952527265082`,0.008212512360652976`,0.006562014679639082`,-0.005211121151164438`}[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"]
0 2 4 6 0 10 20
0 2 4 6 0 10 20
0 2 4 6 0 10 20
ySol[x]
yAp[x]

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

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
\[c(x) y(x)+b(x) (y')[x]+a(x) (y'')[x]=f(x)\]
\(y(0)=0\)
\(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]
\[(\frac{3}{2}) (E^{(\frac{1}{4}) (x^{2})}) y(x)+(\frac{1}{2}) (E^{(\frac{1}{4}) (x^{2})}) x (y')[x]+(E^{(\frac{1}{4}) (x^{2})}) (y'')[x]=(\frac{1}{2}) (E^{(\frac{1}{4}) (x^{2})}) f(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}]
-0.024574622586842638`
0 2 4 6 0 0.5

Στόχος είναι να προσδιορίσουμε την λύση της Σ.Δ.Ε. ως ανάπτυγμα ιδιοσυναρτήσεων. Έτσι, αν $λ_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]
154.7221231584296`

Η προσεγγιστική μας λύση είναι:

yAp[x_] := Sum[yn[j]*eigF[[j]][x], {j, 1, n0}] 
yAp[2]
yn[2]
121.8871842751746`
154.7221231584296`

Θα την αντιπαραβάλουμε με τη λύση του 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"]
0 2 4 6 0 100 200 yAp
0 2 4 6 0 100 200 ySol 
Plot[{ySol[x], yAp[x]}, {x, 0, L}, PlotLegends -> "Expressions"]
0 2 4 6 0 100 200
ySol[x]
yAp[x]