Python-koodia. Huom: alla-oleva ei ole sovelias lapsille.

Turbulenttinen nopeus

import sympy as sym
from sympy import pprint

v = sym.Symbol('v')
F,alpha, beta = sym.symbols('F alpha beta')
a,m,t = sym.symbols('a m t')

eq1L = m*sym.diff( v(t), t)
eq1R = F - alpha*v(t) - beta*v(t)**2
sol1 = sym.dsolve( eq1L - eq1R, v(t) )
pprint( sol1 )

josta tulee jotain. Sympylläkin on pieniä hankaluuksia hahmottaa logaritmin ja eksponenttifunktion laskusääntöjä — symboliset laskimet tuppaavat olemaan liian tarkkoja.

 

Siispä ratkaistaan yllä olevasta separoituvasta yhtälöstä $v(t)$

dv,dt = sym.symbols('dv dt')

eq2L = m* dv(t)/dt
eq2R = F - alpha*v(t) - beta*v(t)**2

sol2 = sym.solve( eq2L - eq2R, dt )

jossa nyt vähän ovelasti ratkaistaan $dt$ infitesimaalin $dv$ avulla. Saadaan kuitenki sama kuin alkuperäisessä:

$dt = \frac{m dv(t)}{F – \alpha v(t) – \beta v^2(t)}$

joka sitten integroidaan puolittain:

$\int_0^tdt = \int_0^v\frac{m dv'(t)}{F – \alpha v'(t) – \beta v’^2(t)}$

jossa tuli pieni muuttujan vaihdos $v \to v’$ (ei ole derivaatta). Vasen puoli on helppo derivoida käsinkin

sym.Integrate(1, (t, 0,t))

Oikean puolen integrointi käsin taas. Suoraan pythonilla. Sijoitetaan ensin $dv$:n paikalle $1$ ja $v(t)$:n paikalle $x$.

x = sym.Symbol('x')
int = sol2[0].subs(dv(t),1).subs(v(t), x)
sym.integrate( int, (x, v,0) )

Saatiin sama hyvin monimutkaisen näköinen ratkaisu kuin edellä dsolvella. Sym.simplify-komento ottaa yhden tekijän yhteiseksi $m$:n kaveriksi. Siitä siis pitää ratkaista $v$. Solve antaa karun vastauksen

NotImplementedError: multiple generators [log(-2*F*sqrt(1/(4*F*beta + alpha**2)) - alpha**2*sqrt(1/(4*F*beta + alpha**2))/(2*beta) + alpha/(2*beta) + v), log(2*F*sqrt(1/(4*F*beta + alpha**2)) + alpha**2*sqrt(1/(4*F*beta + alpha**2))/(2*beta) + alpha/(2*beta) + v)]
No algorithms are implemented to solve equation log(-2*F*sqrt(1/(4*F*beta + alpha**2)) - alpha**2*sqrt(1/(4*F*beta + alpha**2))/(2*beta) + alpha/(2*beta) + v)

 

Positiivisuusoletus

Palataan takaisin alkuun, koska homma ei oikein pelittänyt. Muuttujat ovat positiivisia, joten annetaan Sympynkin se tietää

import sympy as sym
from sympy import pprint

v = sym.Symbol('v', positive=True)
F,alpha, beta = sym.symbols('F alpha beta', positive=True)
a,m,t = sym.symbols('a m t', positive=True)

eq3L = m*sym.diff( v(t), t)
eq3R = F - alpha*v(t) - beta*v(t)**2
sol3 = sym.dsolve( eq3L - eq3R, v(t) )
pprint( sol3 )

Nyt kaikki toimii. Python antaan vastauksen

$m \left(\frac{\log{\left (\operatorname{v}{\left (t \right )} + \frac{- 4 \frac{F \beta}{\sqrt{4 F \beta + \alpha^{2}}} – \frac{\alpha^{2}}{\sqrt{4 F \beta + \alpha^{2}}} + \alpha}{2 \beta} \right )}}{\sqrt{4 F \beta + \alpha^{2}}} – \frac{\log{\left (\operatorname{v}{\left (t \right )} + \frac{4 \frac{F \beta}{\sqrt{4 F \beta + \alpha^{2}}} + \frac{\alpha^{2}}{\sqrt{4 F \beta + \alpha^{2}}} + \alpha}{2 \beta} \right )}}{\sqrt{4 F \beta + \alpha^{2}}}\right) = C_{1} – t$

Tämä jo näyttää vähän $\tanh^{-1}$-funktiolta, mitä sen pitäisi muistuttaa. Nyt pitäisi ratkaista $v(t)$, mutta Sympy ei vieläkään osaa ratkaista yhtälöä. Kokeillaan käsin auttamalla.

Alkuarvo

Käskyt .lhs() ja .rhs() palauttavat kuten halutaan. Siis jaetaan ensin massalla $m$ puolittain ja kerrotaan jakajalla $\sqrt{4F\beta + \alpha^2}$ ja vielä korotetaan eksponenttiin. Saadaan

rhs1 = sym.exp( sym.sqrt(4*F*beta+alpha**2)*sol3.rhs()/m )

Vastaavasti vasemmalle puolelle saadaan järisyttävä tulos

lhs1 = sym.simplify( sym.exp( sym.sqrt(4*F*beta+alpha**2)*sol3.lhs()/m ) )

eli saadaan

$\frac{- 4 F \beta – \alpha^{2} + \alpha \sqrt{4 F \beta + \alpha^{2}} + 2 \beta \sqrt{4 F \beta + \alpha^{2}} \operatorname{v}{\left (t \right )}}{4 F \beta + \alpha^{2} + \alpha \sqrt{4 F \beta + \alpha^{2}} + 2 \beta \sqrt{4 F \beta + \alpha^{2}} \operatorname{v}{\left (t \right )}}$

$=$

$e^{\frac{\left(C_{1} – t\right) \sqrt{4 F \beta + \alpha^{2}}}{m}}$

Nyt, joko saadaan ratkaistuksi $v$.

sol4 = sym.solve( lhs1 - rhs1, v(t)  )

Ratkaisu on vielä pitkä ja hankala, joten määritetään uusia muuttujia ja haetaan yhteisiä tekijöitä

x,y = sym.symbols('x y')
sol5 = sym.collect( sym.simplify( sol4[0] ).subs(4*F*beta, x).subs(sym.sqrt(alpha**2+x)/m,y),(x, alpha) )

$x = 4F\beta$
$y = \frac{\sqrt{\alpha^2 +x}}m$
$v(t) = \frac{\alpha^{2} \left(- e^{C_{1} y} – e^{t y}\right) + \alpha \left(- \sqrt{\alpha^{2} + x} e^{C_{1} y} + \sqrt{\alpha^{2} + x} e^{t y}\right) + x \left(- e^{C_{1} y} – e^{t y}\right)}{2 \beta \sqrt{\alpha^{2} + x} \left(e^{C_{1} y} – e^{t y}\right)}$

Yksinkertaistetaan funktiota määrittämällä alkuehto. Tiedetään, että $v(0) = 0$, josta voidaan ratkaista vakio $C_1$. Sijoitetaan se vielä ja sievennetään. Vastaus tulee heti

C1,C2 = sym.symbols('C1 C2')
C2=sym.solve( sol5.subs(t,0), C1 )
sol6 = ( sol5.subs(C1, C2[0]) )
sol7 = sym.simplify( sol6 ).collect(sym.exp(t*y))

eli vakiolle $C_1$ saadaan

$C_1 = \frac{\log{\left (\frac{- \alpha^{2} + \alpha \sqrt{\alpha^{2} + x} – x}{\alpha^{2} + \alpha \sqrt{\alpha^{2} + x} + x} \right )}}{y}$

ja nopeudelle $v(t)$ tulee ällistyttävän yksinkertainen muoto ottamalla vielä $e^{ty}$ yhteiseksi tekijäksi:
$v(t) = \frac{x \left(e^{t y} -1\right)}{2 \beta \left(\alpha e^{t y} – \alpha + \sqrt{\alpha^{2} + x} e^{t y} + \sqrt{\alpha^{2} + x}\right)}=\frac{x \left(e^{t y} -1\right)}{2 \beta \left(- \alpha + \left(\alpha + \sqrt{\alpha^{2} + x}\right) e^{t y} + \sqrt{\alpha^{2} + x}\right)}
$
Selvitetään vielä, pitääkö alkuperäisessä artikkelissa olleet kertoimet paikkansa. Kertojasta ja jakajan yhteisestä tekijästä saadaan
$AB_p = -\frac{x}{2\beta}$
Jakajalle (ilman tekijää $2\beta$) saadaan

jakaja =sym.fraction( sol7*2*beta )[1]
jajaja.expand().args

ja palautetaan siitä lista, jossa näkyy kaikki kertoimet.
Saadaan lista

$\begin{pmatrix} \sqrt{\alpha^{2} + x}, & \sqrt{\alpha^{2} + x} e^{t y}, & \alpha e^{t y}, & – \alpha \end{pmatrix}$,

josta

$A=\sqrt{\alpha^{2} + x}- \alpha$

$B=\sqrt{\alpha^{2} + x}+ \alpha$

ja saadaan $AB = x$; Melkein sama kuin paperissa.