Rozważmy równanie: $$\frac{dx}{dt} = f(x,t).$$
Nie w chwili $t$ funkcja szukana $x$ będzie miała wartość $x(t)$. Wtedy można rozwinąć ją w szereg Taylora:
$$ x(t+h)=x(t)+h x'(t) +\frac{h^2}{2!} x''(t)+ O(h^3).$$
Pochodne $x'(t),x''(t),...$ można obliczyć korzystając z wyjściowego równania różniczkowego.
Urywając szereg Taylora w powyższym rozważaniu na liniowym członie otrzymujemy schemat pierwszego rzędu:
$$ x(t+h)=x(t)+h f(x(t),t)$$
Obliczmy pochodną równania różniczkowego po czasie, nie zapominając, że $x=x(t)$ jest funkcją czasu:
$$ \frac{d^2}{dt^2}x(t)= \frac{d}{dt}f(t,x(t)) = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}\frac{d x(t)}{dt}= \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}f(t,x(t)) $$
Z drugiej strony rozwijając wyrażenie w szereg Taylora dwóch zmiennych:
$$ f(t+h,x+fh) = f(x,t) + h \frac{\partial f(t,x)}{\partial t}+ h f \frac{\partial f}{\partial x}+ O(h^2)$$
i wstawiając do
$$ x(t+h)=x(t)+h f(t,x(t)) +\frac{h^2}{2!} x''(t)+ O(h^3)=$$
$$=x(t)+h f(t,x(t)) +\frac{h^2}{2}\left[\frac{\partial f}{\partial t} + \frac{\partial f}{\partial x}f(t,x(t))\right]+ O(h^3)=$$
$$=x(t)+\frac{h}{2} f(t,x(t)) +\frac{h}{2}\left[f(t,x) +h\frac{\partial f}{\partial t} + h f \frac{\partial f}{\partial x}\right]+ O(h^3) .$$
otrzymujemy:
$$ x(t+h)=x(t)+\frac{h}{2} f(t,x(t))+\frac{h}{2} f(t+h,x(t)+h f(t,x(t))). $$
Jest to formuła drugiego rzędu, jednak w przeciwieństwie do metody rozwinięcia w szereg Taylora nie wymaga obliczania pochodnej prawej strony równania różniczkowego $f(t,x)$. Wzór ten jest zwany metodą Rungego-Kutty drugiego rzędu.
In [46]:
var('t')
x=function('x',t)
f=function('f')
f(x,t).diff(t)
Out[46]:
In [59]:
Out[59]:
Rozważmy równanie $$ m \dot v = -m g -\gamma v.$$
Dzieląc je przez $mg$ trzymujemy: $$ \frac{1}{g} \frac{dv}{dt} = -1 -\frac{\gamma}{mg}v$$.
Zauważmy, że $$[\frac{\gamma}{mg}]=[\frac{s}{m}],$$ więc dposując kilka członów pod pochodną można napisać:
$$ \frac{d \frac{\gamma}{mg}v}{d \frac{\gamma}{m}t} = -1 - \frac{\gamma}{mg}v.$$
Ponieważ $\frac{\gamma}{m}$ ma wymiar odwrotności sekundy można napisać równanie wyjsciowe w postaci:
$$ \dot u =-1-u,$$ gdzie czas i prędkość są bezwymiarowe i są związane w wielkościami fizycznymi przez:
$$ t'=\frac{\gamma}{m}t$$
$$ u=\frac{\gamma}{mg}v$$
In [58]:
Out[58]:
In [47]:
Out[47]:
Zcałkować równanie $$\dot x = x (1-x) $$ dla $t=(0,4)$ z krokiem $h=0.2$ metodą pierwszego i drugiego rzędu. Porównać wynik z rozwiązaniem analitycznym.
In [1]:
var('x')
f(x)=x*(x-1)
In [2]:
x1=[0.1]
x2=[0.1]
czas=[0.0]
h=0.2
for i in range(40):
x1.append( x1[-1]+f(x1[-1],0)*h )
x2.append( x2[-1]+f(x2[-1],0)*h + 0.5*h^2*f(x2[-1])*f.diff(x)(x2[-1]) )
czas.append(czas[-1]+h)
In [3]:
sol=desolve_odeint(f,0.1,srange(0,10,0.1),x)
p0=line(zip(srange(0,4,0.1),sol[:,0]),figsize=4,color='gray')
p1=point(zip(czas,x1))
p2=point(zip(czas,x2),color='red')
In [4]:
p0+p1+p2
Out[4]:
In [5]:
var('c t')
X = function('X')(t)
desolve(diff(X,t)==f(X), X)
Out[5]:
In [6]:
solve( (X - 1)/X == c*exp(t),X)
Out[6]:
In [7]:
solve( -1/(c*e^0 - 1)==0.1, c)
Out[7]:
In [8]:
plt_exact=plot( -1/(c*e^t - 1).subs(c==-9), (t,0,4),color='green',figsize=4)
In [9]:
p1+p2+plt_exact
Out[9]:
In [57]:
Out[57]:
W systemie Sage istnieje kilka możliwości rozwiązywania równań różniczkowych zwyczajnych (i układów równań pierwszego rzędu).
In [10]:
desolvers?
In [13]:
var('t')
x = function('x', t)
desolve(diff(x,t)==1+x^2, x)
Out[13]:
In [18]:
Out[18]:
In [14]:
x,y=var('x y')
P=desolve_rk4(y^2+1,y,ics=[0,0],ivar=x,output='slope_field',end_points=1.5,step=0.01,thickness=3)
P
Out[14]:
In [23]:
Out[23]:
In [24]:
Out[24]:
In [ ]:
In [28]:
Out[28]:
Metoda stosuje procedurę lsode z odepack (z netlib).
In [26]:
Out[26]:
In [21]:
x,y = var('x,y')
f = [x*(1-y),-y*(1-x)]
tpkt = srange(0,10,0.1)
sol = desolve_odeint(f, [0.5,2],tpkt,[x,y])
p=line(zip(sol[:,0],sol[:,1]),figsize=4)
p
Out[21]:
In [22]:
point(zip(tpkt,sol[:,0]),figsize=4)+\
point(zip(tpkt,sol[:,1]),color='red')
Out[22]:
In [23]:
sol.shape
Out[23]:
In [17]:
p+plot_vector_field(f,(x,0,3),(y,0,2.7))
Out[17]:
In [32]:
y = var('y')
epsilon = 0.001
f = y^2*(1-y)
ic = epsilon
t = srange(0,2/epsilon,1)
sol = desolve_odeint(f,ic,t,y ,rtol=1e-9,atol=1e-10,compute_jac=True)
p = points(zip(t,sol))
p.show(figsize=4)
In [23]:
var('x y')
ratio=200/1.4
plot_vector_field([1,f],(x,0,200),(y,0,1.4),figsize=(6,6),aspect_ratio=ratio) + p
Out[23]:
In [34]:
var('x,y,vx,vy')
f = [vx,vy,0.+0*x,-10.+0*x]
sol = desolve_odeint(f,[0.,0.,3*cos(30/180*pi),3*sin(30/180*pi)],srange(0,.3,0.01),[x,y,vx,vy])
p = point(zip(sol[:,0],sol[:,1]),figsize=4)
p
Out[34]:
In [35]:
var('t')
plt_exact=parametric_plot( ( 3*cos(30/180*pi)*t,3*sin(30/180*pi)*t-10.0*t^2/2 ) , (t,0,2*3*sin(30/180*pi)/10 ) ,color='green')
plt_exact
Out[35]:
In [26]:
(2* 3*sin(30/180*pi)/10)
Out[26]:
In [27]:
cos(30/180*pi)
Out[27]:
In [28]:
tend=.3
N=50.
dt=tend/N
x0,y0,vx0,vy0 = 0.,0.,3*cos(30/180*pi).n(),3*sin(30/180*pi).n()
print 0.5*(vx0^2+vy0^2) + 10*y0
traj=[[x0,y0,vx0,vy0]]
for i in range(N):
x,y,vx,vy=x0+vx0*dt,y0+vy0*dt,vx0,vy0-10.*dt
traj.append([x,y,vx,vy])
x0,y0,vx0,vy0= x,y,vx,vy
plt_euler=point([(tr_[0],tr_[1]) for tr_ in traj],figsize=6,color='red')
print 0.5*(vx^2+vy^2) + 10*y
plt_euler+plt_exact
Out[28]:
In [29]:
# Euler full implicit
tend=.3
dt=tend/N
x0,y0,vx0,vy0 = 0.,0.,3*cos(30/180*pi).n(),3*sin(30/180*pi).n()
traj=[[x0,y0,vx0,vy0]]
for i in range(N):
x,y,vx,vy=x0+vx0*dt, -10*dt^2 + dt*vy0 + y0,vx0,vy0-10.*dt
traj.append([x,y,vx,vy])
x0,y0,vx0,vy0= x,y,vx,vy
plt_euler=point([(tr_[0],tr_[1]) for tr_ in traj],figsize=6,color='red')
print 0.5*(vx^2+vy^2)+y*10
plt_euler+plt_exact
Out[29]:
In [30]:
var('x0,y0,vx0,vy0,x,y,vx,vy,dt')
solve( [x==x0+vx*dt,y==y0+vy*dt,vx==vx0,vy==vy0-10.*dt], [x,y,vx,vy] )
Out[30]:
In [31]:
# Euler Cromer (semi-implicit z wagą 1/2)
tend=.3
dt=tend/N
x0,y0,vx0,vy0 = 0.,0.,3*cos(30/180*pi).n(),3*sin(30/180*pi).n()
print 0.5*(vx0^2+vy0^2)+y0*10
traj=[[x0,y0,vx0,vy0]]
for i in range(N):
vx,vy = vx0,vy0-10.*dt
x,y = x0+0.5*(vx+vx0)*dt,y0+0.5*(vy+vy0)*dt
traj.append([x,y,vx,vy])
x0,y0,vx0,vy0= x,y,vx,vy
plt_euler=point([(tr_[0],tr_[1]) for tr_ in traj],figsize=6,color='red')
print 0.5*(vx^2+vy^2)+y*10
plt_euler+plt_exact
Out[31]:
In [41]:
# Euler Cromer (semi-implicit) oscylator Harmoniczny - zachowanie energii
N=97
omega=0.46
tend=2*pi.n()/omega
dt=tend/N
x0,v0 = 1.1,1.
E0= 0.5*(v0^2+omega^2*x0^2)
traj=[[x0,v0]]
for i in range(N):
v = v0-omega^2*x0*dt
x = x0+v*dt
traj.append([x,v])
x0,v0= x,v
plt_euler=line([(tr_[0],tr_[1]) for tr_ in traj],figsize=6,color='red')
print E0,0.5*(v^2+omega^2*x^2)
var('x,v')
plt_exact=implicit_plot( 0.5*(v^2+omega^2*x^2)==E0,(x,-3,3),(v,-2,2),color='grey' )
plt_euler+plt_exact
Out[41]:
In [43]:
Out[43]:
In [42]:
Out[42]:
In [37]:
Out[37]:
GNU Scientific library – biblioteka funkcji obliczeniowych i naukowych dla C i C++ dostępna na zasadach GPL. Biblioteka jest częścią Projektu GNU.
In [33]:
def f_1(t,Y):
x,y=Y[0],Y[1]
return [x*(1-y),-y*(1-x)]
T = ode_solver()
T.algorithm="rk4"
T.function=f_1
T.ode_solve(y_0=[0.5,2.0],t_span=[0,14],num_points=300)
In [34]:
line( [x_[1] for x_ in T.solution] ,figsize=3)
Out[34]:
$dot v$
In [61]:
Out[61]:
In [60]:
Out[60]:
In [7]:
Out[7]: