Differntialgleichung in c# - Wie?

Hallo Ihr

Und schon das nächste Problem: Ich muss zur Auswertung von Messdaten eine Ableitung ersten Grades hinbekommen (will aus der Änderung Strecke/Zeit eine Geschwindigkeit berechnen). Meine Daten liegen als eindimensionales Array vor.

Wie kann man Differnetialrechnung in C# integrieren? Hab erneut Google bemüht, hier aber nichts brauchbares gefunden...

Gruß
Beren

Ich hab mal als kleines Projekt sowas geschrieben, allerdings weiß ich nicht, ob das dass richtige für dich ist.
Ich habe einfach bestimmte regeln verfasst und die gleichung per regex aufgeteilt und dann die ganzen Regeln darauf angewandt. Ansonsten, wie liegen die Daten vor? Hast du die Funktion daraus schon? Ich mein, wenn du weißt dass es eine Gerade ergibt, kannst du ja eine Gerade bilden und dann einfach die Steigung nehmen (ist die Ableitung).
Also sobald die Funktion bekannt ist, brauchst du nicht die Aufspaltung etc. also das umständliche, sondern kannst direkt spekulieren.

//EDIT: Zum Thema Gleichung aufteilen ist die Postfixschreibweise einer Gleichung hilfreich.

Hab erneut Google bemüht, hier aber nichts brauchbares gefunden... Wenn man weiß nach was man suchen muss ist das anders 🙂Numerische Differentiation

Hi

Nein, es sind Messdaten einer Bewegung. Die Bewegung ist nicht linear. Das heißt, v ist nicht konstant.

Ich habe Daten in einem eindimensionalen Array. Die Daten beschreiben den Abstand zwischen zwei Stauchrahmen einer Schweißmaschine (Stahlindustrie) in mm.

Der Abstand ändert sich während der Schweißung. Ich möchte jetzt die Steigung in jedem Punkt der Abstandskurve haben, um die Geschwindigkeit der Stauchrahmen zu ermitteln.
Leider liegt mein Abi 13 Jahre zurück. Ich weiß zwar noch sehr gut, wofür Differentialrechnung gut ist, bekomme das mathematische aber nicht mehr auf die Reihe.

Gruß
Beren

Jetzt gibts noch ein Problem der Stetigkeit, denn nicht jede Funktion ist auch differenzierbar. Wenn die Daten manchmal Ausreißer haben, musst du zuvor noch eine Regressionsfunktion erstellen, damit sie auch integrierbar ist. Numerik ist ein gutes Stichwort, aber wenn ich ehrlich bin, hatte ich das noch nicht^^. Kommt noch.
Aber es dürften sich auf jeden Fall Klassen finden, die diese Funktionen schon eingebaut haben. Wenn ich was gefunden hab, sag ich Bescheid.

Ich danke Dir. Bin auch gerade kräftig am googlen, hab aber immer noch nichts gescheites gefunden, bzw. ich verstehs nicht (mehr).

Es sind übrigens Messdaten, mit kurzen Ausreissern. Also steckt keine Funktion hinter der Kurve.

Ich brauch einfach nur die Steigung in jedem Punkt. Mal sehen ob ichs selber wieder hinbekomme... Ich schnapp mir jetzt Papier und Stift und versuche mich mal dran.

Gruß
Beren

Mittlerweile wurde ja schon auf den Wikipediaartikel verwiesen, aber vielleicht hilft dir das Folgende trotzdem noch etwas:

Da die Ableitung ein Grenzprozess (lim (||h||->0) (f(t+h)-f(t))/h) ist es nur sinnvoll, davon auf einer nicht-diskreten Menge zu sprechen, was bei deinen endlichen Stützpunkten nicht der Fall ist.

Eine Möglichkeit ist es, durch die Stützstellen zu interpolieren und dann von der interpolierenden Kurve die Ableitung. Welche Form der Interpolation du verwendest hängt aber von den Daten ab. (im Grunde solltest du an Hand von Beispielmesswerten verschiedene Interpolationen ausprobieren und gucken, welche Interpolation am besten dem tatsächlichen Verlauf entspricht)

Als Interpolationsmöglichkeiten gibt es: Polynominterpolation,Splineinterpolation, Interpolation mittels trigonometrischen Polynomen (entspricht diskreter Fouriertransformation, besonders geeignet für periodische Funktionen).
Vermutlich gibt es bereits Implementationen der verschiedenen Interpolationen in C#/.NET. Kannst ja dann etwas experimentieren und gucken, was am besten passt.
Den Wert, den für die Ableitung erhälst, hängt entscheidend davon ab, welche Interpolation du verwendest, du musst also eine verwenden, die das widerspiegelt, was dem tatsächlichen Verlauf entsprechen würde.

Wenn du hingegen nur die Durchschnittgeschwindigkeit in einem Intervall brauchst kannst du einfach die Steigung berechnen also (y1-y0)/(t1-t0)

Der Wikipediaartikel Numerische Differentiation bezieht sich auf Folgendes:

Hast du hingegen eine konkrete Funktion vorliegen, so kannst du die Ableitung mittels kleinem h approximieren durch die Geradensteigung: (f(x+h)-f(x))/h
Allerdings steigt für kleine h der durch Rundungsfehler verursachte Fehler (stell dir ein Epsilon vor, dass im Zähler als Störung dazuaddiert wird, dieses wird durch Teilung mit kleinem h sehr groß)

Oder du berechnest für verschiedene h die Geradensteigung, interpolierst die daraus resultierenden Werte und nimmst den Wert der Funktion bei 0 (das heißt dann Extrapolation).

Edit: Vielleicht könntest du auch ein paar Graphiken von Beispielverläufen online stellen, dann könnte man besser spekulieren, was sich eignet.

Aber Ausreiser sind fürs Ableiten schon mal nicht sehr gut.
Stetigkeit liegt aber auf jedem Fall vor, da es sich ja um eine Bewegung handelt und nicht zwischendrin teleportiert wird 😉

Mal ein paar Verlinkungen:

Algorithm Library .Net
http://www.alglib.net/

CodeProject Solving differntial equations
http://www.codeproject.com/KB/recipes/SolvingDifferentialEqns.aspx

Und Google (für den Rest):
http://www.google.de/search?hl=de&q=differential+c%23&btnG=Google-Suche&meta=

Habs mit Papier und Stift gelöst.

Die Messdaten liegen in einer Auflösung von 1ms vor. Ich mach einfach folgendes:

v = (Abstand [counter+1] - Abstand [counter-1])/0.002

Das ergibt die Steigung in dem Punkt Abstand[counter].

Man... Hätte ich mal früher Blatt und Stift zur Hand genommen....

Gruß
Beren

Dir ist klar, dass v den Abstand zwischen 3 Punkten berechnet? counter+1 und counter-1 und counter 😉
Erklär mir mal bitte die Gleichung, ich versteh nämlich den Sinn davon noch nicht so ganz g

Mhh... schaut doch gut aus, was Beren gemacht hat.

er nimmt als Approximation für den Differentialquotionten an der Stelle x_n den Anstieg zwischen x_n-1 und x_n+1

Man nimmt ja deshalb die umliegenden Werte, weil die Approximation eigentlich für die x-Stelle genau in der Mitte liegt.

Würde er (x_n+1 - x_n) verwenden, wäre das gültig für die Stelle x_n+0,5, aber son Array hat ja eben nur ganzzahlige Indizes 😉

Trotzdem wärs schön, wie benedictbaur schon gesagt hat, wenn du mal ne Beispielskurve posten könntest.
Weil mit ner guten Ausgleichskurve könntest du noch bessere Resultate erzielen.
Schließlich verbietet dir ja keiner Hintergrundwissen deines Prozesses (es kommt oft zu abrupten Änderungen, bei jeder Bewegungsrichtungsänderung kommt es zu einem "Überschwingen" des Schweißarms) mathematisch mit reinzumodellieren, was das Ergebnis natürlich besser macht.

beste Grüße
zommi

PS @ Beren:
Die Differenziation ist natürlich auch ein Filter, und dein Algo ist ein nicht rekursives Filter 😉

PPS: anstatt durch 0.002 zu teilen, isses evtl besser mit 5000 zu multiplizieren.
Wenn ich mich nicht täusche, kann ja eine Multiplikation imernoch schneller ausgeführt werden als ne Division.

Okay, hier das Bild (mein Algo funktiojniert übrigens! 😄)

Gruß
Beren

Hallo,

ich wollte dir kurz mal positives Feedback geben:
Deine Kurve sieht wirklich sehr glatt aus und weist keine erkennbaren Ausreißer auf. Das macht die Sache natürlich einfacher. Du nimmst an, dass alle Messwerte tatsächliche Interpolationspunkte deiner kontinuierlichen Modellfunktion sind. Die Messwerte bilden also eine Diskretisierung deiner Funktion. Dann verwendest du die Finite-Differenzen-Methode bzw. die zentralen Differenzenquotienten zur Approximation der ersten Ableitung. Für deine Gitterweite h=0.001 ist ein Diskretisierungsfehler von O(h^{2) also ca. 10}(-6) zu erwarten.

Das Problem: Sollten wider Erwarten in einigen Messungen doch durch Messfehler hervorgebrachte Ausreißer auftreten, ist die idealisierte Annahme, dass Messpunkte auch Interpolationspunkte sind, hinfällig und der Fehler deiner Methode kann signifikant steigen. In diesem Fall musst du in erster Linie ein Approximationsproblem lösen. Das kann aber schwierig sein, da man evt. keine geeignete Modellfunktion findet. Aber das lässt sich ja noch diskutieren wenn die Notwendigkeit besteht. 😉

Hallo Kartoffel

Vielen Dank für Dein Feedback. Ich denke, daß es im Normalbetrieb keine Ausreißer (Richtungsänderungen) geben wird. Daher bin ich mit meiner Lösung vorerst zufrieden.

Gruß
Beren

Moin zusammen

Jetzt hab ich ein weiteres Problem. Die Steigung in dem Punkt y(n) berechnete ich ja so: my(n)=y(n-1)-y(n+1)/timebase*2

Jetzt muss ich ein Integral über ein anderes Signal bilden. Die Analysenvorschrift sieht so aus (iba Syntax):

Max(Int(XCutValid ([32:8],[Flashing])))/(XLast ([Flashing])-XFirst ([Flashing]))

Die Signale (Schweißstrom = [32:8]) liegen wieder in einem eindimensionalen Array vor. Den Term XCutValid ([32:8],[Flashing]) berechne ich bereits.

// Flow not per surface berechnen 
count = endBooting;
while (!metasensor.Flashing[count])
{
    count++;
}
int startFlashing = count;
while (metasensor.Flashing[count])
{
    count++;
}
 int endFlashing = count;

count = startFlashing;
flow = 0;
while (count < endFlashing)
{
    // Über folgende Zeile ein Integral bilden. Nur wie?
    metasensor.Flow[count] = (iSchweiss.Data[count]);
    if (metasensor.Flow[count] > flow)
    {
        flow = metasensor.Flow[count];
    }
    count++;
}
flow /= (endFlashing * tbase - startFlashing * tbase);

Jetzt fehlt mir nur das Integral über die einzelnen Messpunkte. Ich hab mich schon mit Stift und Papier bemüht, bekomme es aber nicht auf die Reihe...

Gruß
Beren

Ich habe nicht genau verstanden, worüber du ein Integral bilden willst. Aber ich vermute mal Folgendes: Du hast eine ähnliche Diskretisierung einer Funktion wie in deinem Bild zuvor. Darüber willst du nun das Integral bilden.

Ich denke du bekommst eine Approximation an das Integral mit ausreichend Güte, indem du das Integral der zugehörigen linearen Splineinterpolation bestimmst. Lineare Spileinterpolation bedeutet, dass du die Punkte, die dir zwischen zwei Messwerten fehlen auf eine Gerade zwischen den zwei nächsten Messwerten legst. Das Integral besteht aus dem Flächeninhalt unter dieser Funktion und kann durch Addition von Intervallen berechnet werden.

Das geht so:
Seien deine n+1 Messpunkte y0, y1, y2, ..., yn. Die zugehörigen x-Werte ergeben sich aus dem Gitterabstand h (=0.001?). Also x0 = 0, x1 = h, x2 = 2h, ..., xi = i*h, ...

Das Integral von 0 bis n*h ergibt sich durch Aufsummieren der einzelnen Intervalle. Diese wiederrum lassen sich zerlegen in ein Rechteck und ein rechtwinkliges Dreieck.

A_Rechteck(i) = min(yi, y(i+1))*h
A_Dreieck(i) = |yi-y(i+1)|*h/2
A_Intervall(i) = A_Rechteck(i) + A_Dreieck(i)

Das Intergal ist nun einfach die Summe über alle Intervallflächen:

Int(0, ..., k*h) = Summe(i=0, ..., k) A_Intervall(i)

Super Kartoffel!

Du bist spitze. Ich habs geschnallt. Das ich da nicht selber drauf gekommen bin... 🙁

Vielen Dank für Deine ausführliche Erklärung!

Beren

Hier jetzt der vollständige Algo:

for (count = startFlashing; count < endFlashing + 1; count++)
{
    metasensor.Flow[count] = tbase * iSchweiss.Data[count] + (tbase * (iSchweiss[count + 1] - iSchweiss[count])) / 2;
    if (metasensor.Flow[count] > flow)
    {
        flow = metasensor.Flow[count];
    }
}

Gruß
Beren

Da du weder min noch Betrag benutzt hast gehe ich davon aus, dass deine Messdaten immer monoton steigen?

Kartoffel, Du hast recht... Beim berechnen des Dreiecks muss ich den Betrag nehmen... Danke für den Hinweis!

metasensor.Flow[count] = tbase * iSchweiss.Data[count] + Math.Abs(tbase * (iSchweiss[count + 1] - iSchweiss[count])) / 2;

Gruß
Beren

Und beim Berechnen des Rechtecks musst du das Minimum aus yi und yi+1 nehmen! Im Moment benutzt du immer den linken Wert.

Hi ihrs,

ich würd einfach auf Abs, Min, Max verzichten und es über den Mittelwert berechnen:

metasensor.Flow[count] = tbase * (iSchweiss[count + 1] + iSchweiss[count]) / 2;

Darüber hinaus weist du ja bisher jedem x-Wert nur den Flächeninhalt in seinem umgebenden Streifen zu, du musst doch den Flächeninhalt von x=0 bis zu diesem Wert berechnen.

also immer noch temporär speichern (und aufsummieren)

temp = temp + tbase * (iSchweiss[count + 1] + iSchweiss[count]) / 2   //oder auch mit +=
metasensor.Flow[count] = temp;

Aber weil du jetzt nicht symmetrisch zu deinem Wert die Fläche berechnest, is deine Integralfunktion etwas verschoben.
Also Maxima und Minima werden nicht dort sein, wo Nullstellen waren, sondern etwas nach rechts versetzt.

Wenn dich das nich stört, dann lass es so, sonst guck dir nochmal Simpson-Formel an.

beste Grüße
zommi