Navigacija:

math.ba
linearni slučaj
nelinearni slučaj

Uvod (egzistencija i jedinstvenost)

Neka je dat rubni problem $$y''=F(x,y,y'),\:x\in[a,b],\:y(a)=\alpha,\,y(b)=\beta.$$ Funkcija $F$ je neprekidna na skupu $$D=\{(x,y,y'):a\leqslant x\leqslant b,\,-\infty<y<\infty,\,-\infty<y'<\infty\},$$ i neka su i parcijalni izvodi $F_y,\,F_{y'}$ takodje neprekidni na skupu $D.$ Ako vrijedi

  • $F_y(x,y,y')\geqslant 0$ za svaki $(x,y,y')\in D$;

  • postoji konstanta $C,$ takva da je $\left| F_{y'}(x,y,y')\right|\leqslant C,$ za svaki $(x,y,y')\in D,$

tada dati rubni problem ima jedinstveno rješenje.

Primjer Dat je rubni problem $$ y''+e^{-xy}+\sin y'=0,\:x\in[1,2],\:y(1)=y(2)=0.$$ Vrijedi $$F(x,y,y')=-e^{-xy}-\sin y', $$ te kako je $$F_y(x,y,y')=xe^{-xy}>0 \text{ i } \left| F_{y'}(x,y,y')\right|=|-\cos y'|\leqslant 1, $$ ovaj rubni problem ima jedinstveno rješenje.

Ako funkcija $F$ ima oblik $$F(x,y,y')=b(x)y'+c(x)y+f(x), $$ tada je diferencijalna jednačina $$y''=F(x,y,y')$$

linearna. Sada se uslovi pod kojima postoji jedinstveno rješenje mogu oslabiti.

Ako linearni rubni problem $$y''=b(x)y'+c(x)y+f(x),\:x\in[a,b],\:y(a)=\alpha,\,y(b)=\beta,$$ zadovoljava

  • funkcije $b,c$ i $f$ su neprekidne na $[a,b];$
  • $c(x)\geqslant0$ na $[a,b];$

tada rubni problem ima jedinstveno rješenje.

Metoda centralnih razlika

Rješavanje linearnih rubnih problema

Bez umanjenja opštosti nastavićemo analizu na $[0,1]$ umjesto $[a,b].$ Ovo se lako izvrši funkcijom $x\mapsto \frac{x-a}{b-a}.$

Dat je rubni problem \begin{align} -\varepsilon^2 y''+b(x)y'+c(x)y=f(x),\:y(0)=\alpha,\,y(1)=\beta,\:x\in[0,1],\nonumber \label{rubni1} \end{align} i neka je \begin{align} c(x)\geqslant 0,\:x\in[0,1],\nonumber \label{rubni2} \end{align} gdje je $0<\varepsilon\leqslant 1$ i neka su funkcije $b,c$ i $f$ su dovoljno glatke na $[0,1].$
Dati rubni problem ima jedinstveno rješenje.

Napomena 1: Za rješavanje klasičnih linearnih rubnih problema koristiti vrijednost parametra $\varepsilon=1.$ Dobijena šema može se koristiti i za testiranje numeričkog rješavanja odgovarajućeg singularno-perturbacionog rubnog problema.

Ovaj problem rješavaćemo na ekvidistatnoj mreži, čiji je parametar mreže (dužina koraka) $h=\frac{1}{N}$ (u opštem slučaju je $h=\frac{b-a}{n}$), i $N$ je broj podsegmenata segmenta $[0,1],$ dakle tačke mreže dobijamo na sljedeći način $$x_i=ih,\:i=0,1,\ldots,N, \text{ gde je } x_0=0,\:x_N=1. $$

Metoda konačnih razlika sastoji se u diskretizaciji diferencijalne jednačine koristeći tačke mreže $x_i,$ gdje s nepoznate $y_i,\,i=1,\ldots,N-1,$ približne vrijednost tačnog rješenja u tačkama mreža $y(x_i),$ tj. ($y(x_i)\approx y_i $).

Prvi i drugi izvod aproksimiramo centralnim razlikama $$ y'(x)\approx (D^0y)(x):=\frac{y(x+h)-y(x-h)}{2h},$$ i $$y''(x)\approx (D^{+}D^{-}y)(x)=\frac{y(x+h)-2y(x)+y(x-h)}{h^2},$$ gdje su $(D^{+}y)(x):=\frac{y(x+h)-y(x)}{h},\,(D^{-}y)(x):=\frac{y(x)-y(x-h)}{h}.$

Uzimajući $x=x_i,$ dobijamo aproksimacije prvog i drugog izvod u tačkama mreže $$y'(x_i)\approx\frac{y_{i+1}-y_{i-1}}{2h} \text{ i } y''(x_i)\approx\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}.$$

Uvrštajući prethodne aproksimacije dobijamo diferencnu šemu

$$\left(-\frac{\varepsilon^2}{h^2}-\frac{b_i}{2h} \right)y_{i-1} +\left(c_i+\frac{2\varepsilon^2}{h^2} \right)y_i+\left(\frac{b_i}{2h}-\frac{\varepsilon^2}{h^2} \right)y_{i+1}=f_i,\,i=1,2,\ldots,N-1,$$ gdje su $b_i=b(x_i),\,c_i=c(x_i),\,f_i=f(x_i).$

Sada uvrštavajući vrijednosti indeksa $i=1,2,\ldots,N-1$ u diferencnu šemu dobijamo $N-1$ linearnu algebarsku jednačinu sa $N+1$ nepoznatom, te uvrštavajući $y_0=\alpha$ i $y_N=\beta$ na kraju dobijamo $(N-1)\times (N-1)$ sljedeći sistem algebarskih linearnih jednačina

\begin{align} \left( c_1+\frac{2\varepsilon^2}{h^2}\right)y_1+\left(\frac{b_1}{2h}-\frac{\varepsilon^2}{h^2}\right)y_2\hspace{9.5cm}&=f_1-\left(-\frac{\varepsilon^2}{h^2}-\frac{b_1}{2h} \right)\alpha\nonumber \\ \left(-\frac{\varepsilon^2}{h^2}-\frac{b_2}{2h} \right)y_1+\left( c_2+\frac{2\varepsilon^2}{h^2}\right)y_2+\left(\frac{b_2}{2h}-\frac{\varepsilon^2}{h^2}\right)y_3\hspace{5.75cm}&=f_2\nonumber\\ \left(-\frac{\varepsilon^2}{h^2}-\frac{b_3}{2h} \right)y_2+\left( c_3+\frac{2\varepsilon^2}{h^2}\right)y_3+\left(\frac{b_3}{2h}-\frac{\varepsilon^2}{h^2}\right)y_4\hspace{2cm}&=f_3\nonumber\\ \vdots\hspace{3cm}&\nonumber\\ %\vdots&\nonumber\\ \left(-\frac{\varepsilon^2}{h^2}-\frac{b_{N-2}}{2h} \right)y_{N-3}+\left( c_{N-2}+\frac{2\varepsilon^2}{h^2}\right)y_{N-2}+\left(\frac{b_{N-2}}{2h}-\frac{\varepsilon^2}{h^2}\right)y_{N-1}&=f_{N-2}\nonumber\\ \left(-\frac{\varepsilon^2}{h^2}-\frac{b_{N-1}}{2h} \right)y_{N-2}+\left( c_{N-1}+\frac{2\varepsilon^2}{h^2}\right)y_{N-1}&=f_{N-1}-\left(\frac{b_{N-1}}{2h}-\frac{\varepsilon^2}{h^2}\right)\beta.\nonumber \end{align}

Ovo trodijagonalni sistem $$AX=B,$$ gdje su \begin{gather} A= \left(\begin{array}{ccccccccc} c_1+\frac{2\varepsilon^2}{h^2}& \frac{b_1}{2h}-\frac{\varepsilon^2}{h^2}&0&0&\cdots&0&0&0&0\\ -\frac{\varepsilon^2}{h^2}-\frac{b_2}{2h}&c_2+\frac{2\varepsilon^2}{h^2}&\frac{b_2}{2h}-\frac{\varepsilon^2}{h^2}&0&\cdots&0&0&0&0\\ 0&-\frac{\varepsilon^2}{h^2}-\frac{b_3}{2h}&c_3+\frac{2\varepsilon^2}{h^2}&\frac{b_3}{2h}-\frac{\varepsilon^2}{h^2}&\cdots&0&0&0&0\\ &&&&\vdots&&&&\\ 0&0&0&0&\cdots&-\frac{\varepsilon^2}{h^2}-\frac{b_{N-3}}{2h}& c_{N-3}+\frac{2\varepsilon^2}{h^2}&\frac{b_{N-3}}{2h}-\frac{\varepsilon^2}{h^2}&0\\ 0&0&0&0&\cdots&0&-\frac{\varepsilon^2}{h^2}-\frac{b_{N-2}}{2h}& c_{N-2}+\frac{2\varepsilon^2}{h^2}&\frac{b_{N-2}}{2h}-\frac{\varepsilon^2}{h^2}\\ 0&0&0&0&\cdots&0&0&-\frac{\varepsilon^2}{h^2}-\frac{b_{N-1}}{2h}&c_{N-1}+\frac{2\varepsilon^2}{h^2} \end{array} \right)\nonumber\\ X=\left(\begin{array}{c}y_1\\y_2\\\vdots\\y_{N-2}\\y_{N-1}\end{array} \right)\nonumber\\ B=\left(\begin{array}{c}f_1-\left(-\frac{\varepsilon^2}{h^2}-\frac{b_1}{2h} \right)\alpha\\f_2\\\vdots\\f_{N-2}\\f_{N-1}-\left(\frac{b_{N-1}}{2h}-\frac{\varepsilon^2}{h^2}\right)\beta\end{array} \right)\nonumber \end{gather}

Primjer. Riješiti rubni problem $-\varepsilon^2 y''-x y+\sin x\cdot y=-5x\cos 5x+25\varepsilon^2\sin 5x+\sin x\sin 5x,\quad$$y(0)=0,\,y(1)=\sin 5,\,x\in(0,1),$ (tačno rješenje ovog problema je $y(x)=\sin 5x,$ i vrijednost argumenta $\sin 5$ je u radijanima).

In [2]:
import math 
import numpy as np
from matplotlib import pyplot 
from numpy import linalg
from matplotlib import rc
noviparametri = {'figure.figsize': (15, 7), 'axes.grid': True,
             'lines.markersize': 5, 'lines.linewidth': 1.25,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral'}
pyplot.rcParams.update(noviparametri)

ep=1.0
N=12
h=1/N
alpha=0          # rubni uslovi
beta=np.sin(5)   # rubni uslovi

B=np.zeros(N-1)
A=np.zeros((N-1,N-1))

def b(x):
    return -x
def c(x):
    return np.sin(x)
def f(x):
    return -5*x*np.cos(5*x)+25*ep**2*np.sin(5*x)+np.sin(x)*np.sin(5*x)
def tacno(x):
    return np.sin(5*x)

# funkcija f i tačno rješenje y(x)=x^2, treba promijeniti rubni uslov y(1)=1
#def f(x):
#    return -2*ep**2-2*x**2+x**2*np.sin(x)
#def tacno(x):
#    return x**2

B[0]=f(h)-alpha*(-ep**2/h**2-b(h)/(2*h))
A[0][0]=c(h)+2*ep**2/h**2; A[0][1]=b(h)/(2*h)-ep**2/h**2;

B[N-2]=f((N-1)*h)-beta*(b((N-1)*h)/(2*h)-ep**2/h**2)
A[N-2][N-3]=-ep**2/(h**2)-b((N-1)*h)/(2*h)
A[N-2][N-2]=c((N-1)*h)+2*ep**2/h**2

for i in range(1,(N-2)):
    B[i]=f((i+1)*h)
    A[i][i-1]=-ep**2/h**2-b((i+1)*h)/(2*h)
    A[i][i]=c((i+1)*h)+2*ep**2/h**2
    A[i][i+1]=b((i+1)*h)/(2*h)-ep**2/h**2
    
rjesenje=linalg.solve(A,B)
 
rjesenje=np.append(rjesenje,[beta])
rjesenje=np.append([alpha],rjesenje)    

x2=np.linspace(0,1,N+1)
x3=np.linspace(0,1,1000)

pyplot.plot(x2,rjesenje,'ro')
pyplot.plot(x2,rjesenje,'r-.')
pyplot.plot(x3,tacno(x3),'b-')
pyplot.legend(('numeričko','linearni interpolat','tačno'))
#pyplot.savefig('/home/samir/Dropbox/Python/slika1.pdf', dpi=1000)
pyplot.show()
pyplot.close()

print('Greska je E_n=%2.4e'%max(abs(rjesenje-tacno(x2))))
print('Teo. greska je %2.4e'%(1/(N+1)**2))
Greska je E_n=2.1159e-02
Teo. greska je 5.9172e-03

Rješavanje nelinearnih rubnih problema

Idi na: početak linearni

Rješavanje opšteg nelinearnog rubnog problema $$y''=F(x,y,y'),\:x\in[a,b],\:y(a)=\alpha,\:y(b)=\beta, $$ slično je rješavanju linearnog problema. Sistem jednačina, čijim rješavanjem ćemo dobiti aproksimativna rješenje u tačkama mreže, više neće biti linearan, što je najveća teškoća.

Napomena 2: I u ovom slučaju koristeći funkciju $x\mapsto\frac{x-a}{b-a}$ segment $[a,b]$ možemo zamijeniti sa segmentom $[0,1].$

Podijelimo ponovo segment $[a,b]$ na $N$ jednakih podsegmenata tačkama $x_i=a+ih,\:i=0,1,\ldots,N,\:h=\frac{b-a}{N}.$ Prvi i drugi izvod aproksimiraćemo na isti način kao i u linearnom slučaju $$y'(x_i)\approx\frac{y_{i+1}-y_{i-1}}{2h},\quad y''(x_i)\approx\frac{y_{i+1}-2y_i+y_{i-1}}{h^2},$$ gdje je $y_i$ aproksimacija tačne vrijednosti $y(x_i)$ u tačkama mreže $x_i,\,i=0,1,\ldots,N.$

Koristeći aproksimacije za prvi i drugi izvod dobijamo diferencnu šemu $$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=F\left(x_i,y_i,\frac{y_{i+1}-y_{i-1}}{2h} \right),\,i=1,2,\ldots,N-1.$$

Ponovo kao i linearnom slučaju uvrštavajući vrijednosti indeksa $i$ u diferencnu šemu i uvrštavajući $y_0=\alpha,\,y_N=\beta,$ dobijamo

\begin{align} \frac{\alpha-2y_1+y_{2}}{h^2}&=F\left(x_1,y_1,\frac{y_2-\alpha}{2h} \right)\nonumber\\ \frac{y_1-2y_2+y_{3}}{h^2}&=F\left(x_2,y_2,\frac{y_3-y_1}{2h} \right)\nonumber\\ \quad\vdots\nonumber\\ \frac{y_{N-2}-2y_{N-1}+\beta}{h^2}&=F\left(x_{N-1},y_{N-1},\frac{\beta -y_{N-2}}{2h} \right). \nonumber \end{align} Rješavanjem prethodnog sistema dobijamo aproksimativne vrijednosti $y_i,\,i=1,2,\ldots,N-1,$

Primjer. Izračunati numeričko rješenje rubnog problema $y''=y+y^2-6e^{-x}\cos 3x-9e^{-x}\sin 3x -e^{-2x}\sin^2 3x,\quad x\in[0,10],$$\:y(0)=0,\,y(10)=e^{-10}\sin 30.$

In [15]:
from scipy.optimize import fsolve
import numpy as np
import matplotlib 
from matplotlib import pyplot

from matplotlib import rc
noviparametri = {'figure.figsize': (12, 8), 'axes.grid': True,
             'lines.markersize': 5, 'lines.linewidth': 1.25,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral'}
pyplot.rcParams.update(noviparametri)

n=64; #broj podsegmenata 
alpha=0; beta=np.exp(-10)*np.sin(30);
xa=0;xb=10
pocetni=.5
h=(xb-xa)/n
  
#def F(x,y):
#    F=(-y+y*y-2+x**2-x**4) 
#    return F
#def tacno(x):
#    return x**2

def F(x,y):
    return y+y**2-6*np.exp(-x)*np.cos(3*x)-9*np.exp(-x)*np.sin(3*x)-np.exp(-2*x)*(np.sin(3*x))**2
def tacno(x):
    return np.exp(-x)*np.sin(3*x)

poc=np.zeros(n-1)
for i in range(n-1):
    poc[i]=pocetni
    
x=np.zeros(n+1)
for i in range(n+1):
    x[i]=xa+i*h

def nelin(y):
    shema=np.zeros(n-1)
    shema[0]=alpha-2*y[0]+y[1]-h*h*F(x[1],y[0])
    shema[n-2]=y[n-3]-2*y[n-2]+beta-h*h*F(x[n-1],y[n-2])
    for i in range(1,n-2):
        shema[i]=y[i-1]-2*y[i]+y[i+1]-h*h*F(x[i+1],y[i])
    return shema

rjesenje=fsolve(nelin,poc)
rjesenje=np.append([alpha],rjesenje)
rjesenje=np.append(rjesenje,[beta])

x1=np.linspace(xa,xb,1000)

pyplot.plot(x1,tacno(x1),'b-.-.',x,rjesenje,'ro',x,rjesenje,'ro-.')
pyplot.legend(('Tačno ','Numeričko ','Linearni interpolat'))
pyplot.show()
pyplot.close()

greska=max(abs(rjesenje-tacno(x)))
print('Greska je En=%2.5e'%greska)
print('Teorijska vrijednost greske je Et=%2.5e'%h**2)
Greska je En=9.39907e-03
Teorijska vrijednost greske je Et=2.44141e-02

Literatura

  1. R.L.Burden, J.D. Faires, Numerical analisys, BROOKS/COLE, Pacific Grove, CA, USA, 2001.
  2. H. B. Keller, Numerical methods for two-points boundary value problems, Dover publications, Inc., Mineola, New York, USA, 1992.

Vrati se na početak