[C#] Differentialgleichungen mit Einschrittverfahren lösen

    • C#
    • .NET (FX) 4.5–4.8

    Es gibt 3 Antworten in diesem Thema. Der letzte Beitrag () ist von jvbsl.

      [C#] Differentialgleichungen mit Einschrittverfahren lösen

      -1-
      Hallo,

      mit den nachfolgenden Klassen ist es möglich, Differentialgleichungen zu lösen.
      Dabei sind die Methoden ( RK4, RK5, Euler) beliebig austauschbar:


      Die IApproximation-Klasse

      C#-Quellcode

      1. public abstract class IApproximation
      2. {
      3. public delegate Vector2 f_delegate(double x, Vector2 y);
      4. public f_delegate f { get; set; }
      5. public double Delta { get; set; }
      6. public double Time { get; set; }
      7. public Vector2 Output { get; set; }
      8. public abstract Vector2 Run(Vector2 y);
      9. }


      Die einzelnen Verfahren

      Euler

      C#-Quellcode

      1. public Euler(Vector2 _0, double delta)
      2. {
      3. Output = _0;
      4. Delta = delta;
      5. }
      6. public override Vector2 Run(Vector2 y)
      7. {
      8. Output += Delta * f(Time, y);
      9. Time += Delta;
      10. return Output;
      11. }
      12. }


      RK4

      C#-Quellcode

      1. public class RK4 : IApproximation
      2. {
      3. public RK4(Vector2 _0, double delta)
      4. {
      5. Output = _0;
      6. Delta = delta;
      7. }
      8. public override Vector2 Run(Vector2 y)
      9. {
      10. Vector2 k1 = Delta * f(Time, y);
      11. Vector2 k2 = Delta * f(Time + .5 * Delta, y + .5 * k1);
      12. Vector2 k3 = Delta * f(Time + .5 * Delta, y + .5 * k2);
      13. Vector2 k4 = Delta * f(Time + Delta, y + k3);
      14. Time += Delta;
      15. Output += 1.0 / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
      16. return Output;
      17. }
      18. }


      und RK5

      C#-Quellcode

      1. public class RK5 : IApproximation
      2. {
      3. public RK5(Vector2 _0, double delta)
      4. {
      5. Output = _0;
      6. Delta = delta;
      7. }
      8. public override Vector2 Run(Vector2 y)
      9. {
      10. Vector2 k1 = Delta * f(Time, y);
      11. Vector2 k2 = Delta * f(Time + ((3.0 * Delta) / 5.0), y + ((1 * 3.0) / 5.0) * k1);
      12. Vector2 k3 = Delta * f(Time + ((2.0 * Delta) / 5.0), y + ((4.0 * 1) / 15.0) * k1 + ((2.0 * 1) / 15.0) * k2);
      13. Vector2 k4 = Delta * f(Time + (Delta / 5.0), y + ((3.0 * 1) / 20.0) * k1 + (1 / 20.0) * k3);
      14. Vector2 k5 = Delta * f(Time + ((4.0 * Delta) / 5.0), y + ((-1 / 5.0)) * k1 + ((-2.0 * 1) / 5.0) * k2 + ((7.0 * 1) / 5.0) * k3);
      15. Vector2 k6 = Delta * f(Time + Delta, y + ((59.0 * 1) / 84.0) * k1 + ((40.0 * 1) / 21.0) * k2 + ((-165.0 * 1) / 28.0) * k3 + ((20.0 * 1) / 7.0) * k4 + ((10.0 * 1) / 7.0) * k5);
      16. Output += (k1 / 12.0) + ((25.0 * k3) / 72.0) + ((25.0 * k4) / 144.0) + ((25.0 * k5) / 72.0) + ((7.0 * k6) / 144.0);
      17. Time += Delta;
      18. return Output;
      19. }
      20. }


      Die Vektor2 Klasse:

      C#-Quellcode

      1. public struct Vector2
      2. {
      3. public double X { get; set; }
      4. public double Y { get; set; }
      5. public double Length()
      6. {
      7. return Math.Sqrt(X * X + Y * Y);
      8. }
      9. public Vector2(double x, double y)
      10. {
      11. X = x;
      12. Y = y;
      13. }
      14. public static Vector2 operator *(int a, Vector2 b)
      15. {
      16. return new Vector2(b.X * a, b.Y * a);
      17. }
      18. public static Vector2 operator *(double a, Vector2 b)
      19. {
      20. return new Vector2(b.X * a, b.Y * a);
      21. }
      22. public static Vector2 operator *(Vector2 b, double a)
      23. {
      24. return new Vector2(b.X * a, b.Y * a);
      25. }
      26. public static Vector2 operator /(Vector2 b, double a)
      27. {
      28. return new Vector2(b.X / a, b.Y / a);
      29. }
      30. public static Vector2 operator +(Vector2 a, Vector2 b)
      31. {
      32. return new Vector2(a.X + b.X, a.Y + b.Y);
      33. }
      34. public static Vector2 operator -(Vector2 a, Vector2 b)
      35. {
      36. return new Vector2(b.X - a.X, b.Y - a.Y);
      37. }
      38. public static Vector2 operator -( Vector2 b)
      39. {
      40. return new Vector2(-b.X, -b.Y);
      41. }
      42. public override string ToString()
      43. {
      44. return "X: " + X + ";" + "Y: " + Y;
      45. }
      46. public PointF ToPointF()
      47. {
      48. return new PointF((float)X, (float)Y);
      49. }
      50. }



      Beispiel:
      Integriert wird die Funktion
      y'=y
      mit der Bedingung y(0)=1

      C#-Quellcode

      1. static void Main(string[] args)
      2. {
      3. IApproximation approx = new RK4(new Vector2(0, 1), 0.000001);
      4. approx.f = (x, y) =>
      5. {
      6. return y;
      7. };
      8. while (approx.Time <= 1.0)
      9. approx.Run(approx.Output);
      10. Console.WriteLine(approx.Output);
      11. Console.Read();
      12. }


      Ergebnis:


      Und das ist offensichtlich richtig:
      Denn
      f(x)= e^x
      f'(x)= e^x
      also f'(x)=f(x) oder y'=y

      Dabei ist while(approx.Time <= n) dasselbe wie f(n)=e^n

      und y'=n*y = f(x)=e^(n*x)

      Viel Spaß.

      EDIT: _ Gelernt mit: web.archive.org/web/2017042106….us/post/runge-kutta-in-c
      _
      Und Gott alleine weiß alles am allerbesten und besser.

      Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „φConst“ ()

      Ist float.
      Wieso immer f Suffix am Ende setzen, wenn man double nutzen kann? :)
      Ist halt meine Art.

      Und Gott alleine weiß alles am allerbesten und besser.

      Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „φConst“ ()