• Register

dvojno integriranje

+5 votes
473 views
Zanima me kako se lahko dvojno integrira funkcijo z eno spremenljivko s pomočjo scipy.integrate.dlbquad?

Natančneje, imam izraz za pospešek, neznanka je čas(t), integrirati pa želim od 2 do 0 da dobim izraz za pot.
asked Dec 20, 2015 by tomi (1,180 points)
edited Dec 20, 2015 by tomi

2 Answers

+5 votes
 
Best answer

scipy.integrate.dlbquad je namenjen funkcijam z dvema spremenljivkama f(x, y), torej ni primerna za vaš namen. 

V vašem primeru gre v bistvu za reševanje nedoločenega integrala:

Imate dve možnosti reševanja

  1. Obravnavamo integral F(t) kot diferencialno enačbo: 

      dF/dt = f(t)

      2. Integriranje F(t) čez celotno območje s spreminjajočo zgornjo mejo.

Obe rešitve si poglejmo na primeru:

xdd = sin(t) # pospesek, F(t)

Izračunati želimo hitrosti in pomike od časa 0 do 10 s. To pomeni, da moramo integrirati enkrat oz. dvakrat. Najprej definirajmo začetne pogoje:

x0 = 0  # začetni pomik
xd0 = -1  # začetna hitrost

1. Poglejmo si reševanje z diferencialnimi enačbami:

a) definirajmo vektor stanj (diferencialno enačbo 2. reda pretvorimo na dve diferencialne enačbe prvega reda):

def resevanje_diferencialne_enacbe(state,t):
    # razpkirajmo vektor stanj
    x = state[0]  # pomiki
    xd = state[1]  # hitrosti, 1. differencialna. enačba

    # izracunajmo pospeske
    xdd = np.sin(t)  # 2. differencialna enačba
    # vrnimo resitev 
    return [xd, xdd]
state0 = [x0, xd0]  # Začetni pogoji (pozor: 
# napačni začetni pogoji lahko močno vplivajo na rešitev)
time = np.arange(0.0, 10, 0.001)  # Definirana časovna os
# reševanje
state = integrate.odeint(resevanje_diferencialne_enacbe, state0, time)  
# Izris
plt.plot(time1, state)
plt.xlabel('cas (s)')
plt.ylabel('stanja')
plt.title('Resevanje diferencialne enacbe')
plt.ylim([-1, 1])
plt.legend(('$x$ (m)', '$\dot{x}$ (m/s)'))
 

Vidimo, da je rešitev pravilna!

2.) Sedaj si poglejmo še reševanje še z spreminjajočo se zgornjo mejo:

def func(x):
    return np.sin(x)
t = np.linspace(0, 10, 1001)
n = len(t)
hitrost = np.zeros(n)
pomik = np.zeros(n)
# Določitev hitrosti in pomikov
for i in range(n):  
    hitrost[i] = xd0 + integrate.quad(func, 0, t[i])[0]
    # Gaussove kvadrature ne moremo več uporabiti, saj nimamo podane funkcije
    pomik[i] = x0 + integrate.trapz(hitrost[:i+1], t[:i+1])   
# Izris
plt.plot(t, pomik)
plt.plot(t, hitrost)
plt.xlabel('cas (s)')
plt.ylabel('stanja')
plt.title('Resevanje diferencialne enacbe')
plt.ylim([-1, 1])
plt.legend(('$x$ (m)', '$\dot{x}$ (m/s)'))
 

Vidimo, da je tudi v tem primeru rešitev pravilna.

answered Dec 21, 2015 by blaz (41,550 points)
selected Dec 21, 2015 by blaz
Kakor sem razumel moramo prvič integrirati z nedoločenim integralom, kar pa s scipy.integrate.quad ni mogoče. Šele pri drugem integriranju določimo meje. S katero funkcijo je torej to možno?
Sem popravil odgovor.
+1 vote
Funkcija scipy.integrate.dlbquad je namenjena za funkcije dveh sprejemljivk (x, y), ki jih integriramo po x in y, zato zahteva še y = gfun(x) in y =hfun(x), ki sta funkciji spodnje in zgornje meje krivulje.

Če imaš pospešek, ki je funkcija ene sprejemljivke rajši uporabi sympy funkcijo integrate(). Prvič narediš nedoločeni integral, da dobiš izraz za hitrost, drugič pa določeni integral v mejah (t, (0,2)) da izračunaš pot.
answered Dec 20, 2015 by zanpirc (1,200 points)
naloga zahteva numerični izračun, sympy pa izračuna samo simbolično, zato ne vem če se to sme uporabiti pri tej nalogi
Simbolično računa integral več kot 1 uro in ga ne uspe izračunati, kar pomeni možen je samo numeričen izračun. Torej z integrate() si ne moremo veliko pomagati
...