Mathe und mehr…

Moin,

hier die fertigen Dateien:

Klick mich!

Gruß, Daniel

from pylab import *
from scipy.optimize import leastsq

fr = arange(0.05,2.05,0.05)
dae = arange(0.2,0.6,0.8)
import pendel_ode

k=0
ampl=zeros(40)

for y in dae:
for x in fr:
pendel_ode.f_d=x
pendel_ode.r=y
dt= 0.005
t= arange(0,20,dt)
start=[20,0]
out= rk4(pendel_ode.calc_p_d,start,t)
[pxx,freq]=psd(out[2000:,0], NFFT=40000, Fs=2000)
ampl[k]=max(pxx)
k=k+1
print(x)
figure()
# Fit starten
fitfunc= lambda w,d: 10/(sqrt((1-w**2)**2+(2*d*w)**2))
errfunc = lambda p,w,y: fitfunc(w,p[0])-y
p, success = leastsq(errfunc, 1 , args=(fr,ampl))
print(success)
plot(fr,fitfunc(fr,p[0]))
plot(fr,ampl)

Moin.

Hier nochmal die Befehlskette für Freitag (Sinus mit 50 Hz, Spektrum erstellen)

# 1000 t-Werte erzeugen, an denen der Sinus berechnet wird.
# Zwischen 0 und 1 ist praktisch,
# dann kann man die 50 Perioden nachzählen.

t = linspace(0,1,1000)

# Wir rechnen einen Sinus mit 50 Hz an den t-Werten au a
# (wir wissen: Sin mit 50 Hz: sin(2*pi*omega*t)

y= sin(2*pi*50*t)

# Das ganze erstmal ansehen

plot(t,y)

# jetzt das Spektrum erzeugen. Der interessante
# Wert ist Fs, die Samplingrate.
# Mit dem sollte man mal spielen! Bei mir gab es einen Peak bei 50 Hz
# wenn ich den Wert auf 1000 setze. NFFT ist (glaube ich)
# nicht so wichtig, sondern ein Parameter für die Fourier-Trafo.

psd(y, NFFT=256, Fs= 1000)

Moin Moin

Ihr bekommt die Quelltexte, die ich heute auch erläutert habe, alle hier:
Klick mich!

Ich werde das Ganze, sobald ich ein wenig Zeit finde, (mit Glück das kommende Wochenende), mit Kommtaren versehen, damit man alles nochmal nachvollziehen kann.

Viele Grüße,

Daniel