Messwerte integrieren mit Simpson-Regel

Mercator

Hi,
Ich bekomme alle paar Millisekunden einen neuen Messwert, den ich gerne aufintegrieren möchte. Dabei wird immer das aktuelle Ergebnis gespeichert und nicht jeder einzelne Messwert, weil es theoretisch unendlich viele sind.

Ich weiß, dass es Integrationsformeln verschiedener Ordnung gibt (Rechteck, Trapez, Simpson, Simpson 3/8, Boderegel, usw.). Meine Probleme fangen ab der 3. Ordnung an (Simpson).
Ich kann mir vorstellen, dass man jeweils 3 Messwerte abwarten kann, dann integriert, 3 neue Messwerte abwartet und wieder integriert.
Wie man allerdings bei jedem neuen Messwert ein aktuelles Integrationsergebnis angeben kann, ist mir nicht klar, da sich die Intervalle - bestehend aus den letzten 3 Messungen - überschneiden müssten und das Ergebnis zu groß ausfallen würde.

In meinem Lehrbuch steht aber
X_k+1 = X_k + ( V_k-1 + 4 V_k + V_k+1 ) * (1/3) dt

Wenn ich nun für jeden neuen Wert die letzten drei Messwerte für die V's einsetze, dann kommt auch tatsächlich etwas zu großes raus.

Was mache ich falsch?

snOOfy

Das wird wohl mit der Simpson-Regel nicht so funktionieren, da man immer über das ganze Intervall integriert.
Ich würde stattdessen vorschlagen: Speichere immer die letzten 4 Punkte (inklusive des aktuellen), lege ein Polynom dritten Grades durch diese 4 Punkte und integriere es jeweils nur im Bereich zwischen den letzten beiden Punkten. Bei einer einfachen Funktion wie einem Polynom kannst du diese Integration analytisch durchführen.

knivil

http://de.wikipedia.org/wiki/Simpsonregel -> Summierte Simpsonregel -> alternative Formulierung
http://de.wikipedia.org/wiki/Runge-Kutta-Verfahren -> Allgemeine Formulierung -> Beispiel

Btw: Bei deiner Formel dividierst du durch 3, ist nicht 6 richtig?

krümelkacker

Als jemand der sich mit digitaler Signalverarbeitung auskennt, würde ich Dir ja gerne helfen, aber ich weiß nicht so recht, was Du eigentlich machen willst. Wenn Du Dich quasi nur für das Endergebnis interessierst, brauchst Du nur all die Abtastwerte aufaddieren. Kompliziertere Interpolationen helfen da nicht wirklich. Wenn Du aus einem Signal ein neues Signal erzeugen willst (dessen Integral), dann benötigst Du einen "Integrator". Ein "Integrator" ist nichts anderes als ein spezielles, digitales Filter; denn Integration und Differentiation lassen sich als Faltung darstellen. Wie dieser Filter aussehen soll, hängt von Deinen Ansprüchen ab. Worauf kommt es Dir da an?

Gibt es denn einen relevanten Punkt warum du den billigsten aller Integratoren nicht verwenden willst?

Mercator

knivil schrieb:

Btw: Bei deiner Formel dividierst du durch 3, ist nicht 6 richtig?

Kommt vermutlich das gleiche raus. Habe die angegebene Formel auch an Beispielen überprüft. Sie ist korrekt.

Tim schrieb:

Gibt es denn einen relevanten Punkt warum du den billigsten aller Integratoren nicht verwenden willst?

Ich integriere Beschleunigungen zweimal zur Position. Eine Interpolation höherer Ordnung sollte nach meinem Verständnis das Bewegungsverhalten besser darstellen. Außerdem bestätigt die Literatur, dass es wenigstens einen kleinen Genauigkeitsvorteil gibt. Höhere Ordnungen des Interpolationspolynoms werden zum Beispiel in Flugzeugen verwendet.

krümelkacker schrieb:

Wenn Du aus einem Signal ein neues Signal erzeugen willst (dessen Integral), dann benötigst Du einen "Integrator". Ein "Integrator" ist nichts anderes als ein spezielles, digitales Filter; denn Integration und Differentiation lassen sich als Faltung darstellen. Wie dieser Filter aussehen soll, hängt von Deinen Ansprüchen ab. Worauf kommt es Dir da an?

Wie gesagt, es handelt sich um Beschleunigungen. Mit einem Interpolationspolynom höherer Ordnung möchte ich einen möglichst natürlichen Verlauf zwischen den Messungen erreichen, um die Genauigkeit der Positionsberechnung zu steigern.
Ein Integrator stellt sich meines Wissens im Frequenzbereich als 1/s dar. Um diesen Filter im zeitdiskreten zu implementieren, müsste man aber wieder eine Approximation durchführen. Mein Regelungstechnik-Skript schlägt zum Beispiel die Tustin-Approximation vor. Für ein Integral scheint mit diese Approximation aber nichts anderes als die Trapezregel zu sein.

snOOfy schrieb:

Das wird wohl mit der Simpson-Regel nicht so funktionieren, da man immer über das ganze Intervall integriert.
Ich würde stattdessen vorschlagen: Speichere immer die letzten 4 Punkte (inklusive des aktuellen), lege ein Polynom dritten Grades durch diese 4 Punkte und integriere es jeweils nur im Bereich zwischen den letzten beiden Punkten.

Das sollte funktionieren. Scheint wohl auch die einzige Möglichkeit zu sein.

snOOfy schrieb:

Ich würde stattdessen vorschlagen: Speichere immer die letzten 4 Punkte (inklusive des aktuellen), lege ein Polynom dritten Grades durch diese 4 Punkte und integriere es jeweils nur im Bereich zwischen den letzten beiden Punkten. Bei einer einfachen Funktion wie einem Polynom kannst du diese Integration analytisch durchführen.

Das ist ja gerade die Simpson-Regel für nur 4 Punkte. Und dann addierst du die Flächeninhalte zusammen, statt über das gesamte Intervall zu integrieren.

Die Frage die sich mir stellt: Wenn du tatsächlich so viele Werte aufaddieren willst, dass du sie nicht mehr Speichern kannst, dann kriegst du ziemlich schnell Probleme mit Rundungsfehlern. Da dürfte dann am Ende nichts rauskommen, was noch mit dem eigentlichen Ergebnis zu tun hat. Eventuell muss man sich da etwas geschickter anstellen, als einfach alle Werte in der Reihenfolge zu addieren, wie sie auftauchen.

knivil

Mercator schrieb:

Kommt vermutlich das gleiche raus. Habe die angegebene Formel auch an Beispielen überprüft. Sie ist korrekt.

Wahrscheinlich liegt dein Problem ganz woanders.

Scheint wohl auch die einzige Möglichkeit zu sein.

Nee, ist es nicht. Warum sollte sukkzessive Integration nicht funktionieren? Kuckst du hier: http://www.youtube.com/watch?v=f594te-zIXM . Es wird auch gesagt, dass ein "Polynom dritten Grades" keine Verbesserung bringt.

simpson.h:

#ifndef SIMPSON_H
#define SIMPSON_H
class Simpson
{
  public:
    Simpson(double h, double a);
    double add(double p1, double p2);
    double curr() const;

  private:
    double h;
    double last;
    double area;
};
#endif

simpson.cpp:

#include "simpson.h"

Simpson::Simpson(double h, double a) : h(h), last(a), area(0)
{
}

double Simpson::add(double p1, double p2)
{
  area += last + 4*p1 + p2;
  last = p2;
  return h*area/3;
}

double Simpson::curr() const
{
  return h*area/3;
}

main.cpp:

#include "simpson.h"
#include <cmath>
#include <iostream>

int main()
{
  double h = 0.1;
  double v = 0;
  Simpson s(h, 0);
  while(v < 100)
  {
    double sample1 = std::sin(v + h);
    double sample2 = std::sin(v + 2*h);

    s.add(sample1, sample2);

    std::cout << s.curr() << " " << 1 - std::cos(v+2*h) << "\n";
    v += 2*h;
  }
  return 0;
}