Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
104 views
Tweet
Share
Share

Čo vieme o danom systéme dif. rovníc (ODE)

použitý free open source CAS softvér - Sage

PC verzia softvéru - http://www.sagemath.org/index.html Cloud verzia softvéru - https://cocalc.com/

autor: Jozef Hanč P.J. Šafárik University in Košice, Slovakia [email protected]

Konkrétny tvar

závislé premenné sú y,ty,t; nezávislá xx; k,yfk, y_f sú parametre:

1.DR:y⋅z(x)⋅k−1=dydx2.DR:−12y=dtdx\Large {\rm 1.DR:}\quad y\cdot z(x)\cdot k-1 =\dfrac{dy}{dx} \\ \Large {\rm 2.DR:}\quad\hspace{1 cm} -\dfrac{1}{2\sqrt{y}} =\dfrac{dt}{dx}

 

kde

  • z(x)=1k⋅yfe−kx\large z(x) = \dfrac{1}{\mathit{k\cdot y_f}}e^{-k x}

  • počiatočné podmienky: y(1)=0,t(1)=0\large y(1)=0, t(1)=0

  • predpoklady o parametroch: k,yf∈R+\large k, y_f \in R^+

Analytické riešenie pre danú počiatočnú podmienku

Vieme zapísať pomocou

  • metódy variácie konštánt pre prvú rovnicu

  • separácie premenných pre druhú rovnicu

v tvare:

y(x)=e−z(x)k[Ei(z(x))−C1]\Large y(x) = \dfrac{e^{-z(x)}}{k}{\left[{\rm Ei}\left(z(x)\right) - C_1\right]}

t(x)=12k∫z1z(x)eu/2uEi(u)−C1du\Large t(x) = \dfrac{1}{2 \sqrt{k}} \displaystyle \int_{z_1}^{z(x)} \dfrac{e^{u/2}}{u\sqrt{{\rm Ei}(u)-C_1}} du

 

kde

  • z1=z(1),C1=Ei(z(1))\Large z_1=z(1), C_1 = {\rm Ei}(z(1))

  • Ei\rm Ei je exponenciálny integrál: Ei(x)=V.P.∫−∞xettdt,x∈R\quad \large{\rm Ei}(x) = {\rm V.P.} \displaystyle\int_{-\infty}^x \dfrac{e^t}{t}dt, \quad x \in R

  • V.P. = Cauchyho hlavná hodnota integrálu

Numerické riešenie aj so smerovým poľom

Z našich dát: k≈5,yf≈1/260k \approx 5, y_f \approx 1/260

Použitie solvera odeint zo SciPy https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html

Odeint má v sebe solver lsoda z FORTRAN library odepack http://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html

# premenné a parametre
%var y,t
yf = 1/260; k = 5
# množina bodov x
dx = 0.001
xh = [1,1-dx..0]
# počiatočná podmienka
ics = [1e-10,0]
# numerické riešenie ODE
F = [y/yf*exp(-k*x)-1,-1/(2*sqrt(y))]
sol = desolve_odeint(F,ics,xh,[y,t],ivar=x)
yh = sol[:,0]; th = sol[:,1]
yvsx = zip(xh,yh)
tvsx = zip(xh,th)

Smerové pole pre 1. DR

# pomocné funkcie
f1(x,y) = ( 1, y/yf*exp(-k*x)-1 )
cp(x) = yf*exp(k*x)
# Rozmery smerového poľa
xmin = 0; xmax = 1
ymin = 0; ymax = 1
pocet = 30
# Smerové pole a presné riešenie
gs = plot_vector_field(-f1(x,y)/norm(f1(x,y)), (x,xmin,xmax), (y,ymin,ymax), plot_points=pocet, color='green')
# krivka numerického riešenia
gy = line(yvsx, color = 'red', legend_label='solution: $y = y(x), y(1)=0 $')
# krivka kritických bodov dy/dx = 0
gc = plot(cp(x), (x, xmin, xmax), color ='blue', legend_label='critical points: $dy/dx = 0$')
# výsledný graf
g = gs+gy+gc
g.axes_labels(['$x$','$y$'])
g.show(figsize=(8,8))

Smerové pole pre 2. DR

# parametre a funkcie
z(x) = 1/(k*yf)*exp(-k*x)
z1 = z(1)
C1 = Ei(z1)
yp(x) = exp(-z(x))/k*(Ei(z(x))-C1)  #analytické riešenie
f2(x,y) = ( 1, -1/(2*sqrt(yp(x))) )
# grafické parametre
xmin = 0; xmax = 1
ymin = 0; ymax = 3.1
pocet = 30
# smerové pole
gs = plot_vector_field(-f2(x,y)/norm(f2(x,y)), (x,xmin,xmax), (y,ymin,ymax), plot_points=pocet, color='green')
# krivka numerického riešenia
gt = line(tvsx, color='darkviolet', legend_label='solution: $t = t(x), t(1) = 0$')
# výsledný graf
g = gs+gt
g.axes_labels(['$x$','$t$'])
g.show(figsize=(8,8))