math.ba
linearni slučaj
nelinearni slučaj
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
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
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.
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}
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 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.
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,$
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