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
- 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.