Реализация метода Рунге-Кутты для параллельного решения системы ОДУ - C#

Узнай цену своей работы

Формулировка задачи:

Доброго времени суток. Коллеги, не найдется ли у кого реализации метода Рунге-Кутты для решения систем ОДУ на C#, который можно было бы распараллелить? Ну или просто реализации, а над параллельностью сам подумаю. Заранее спасибо!)

Решение задачи: «Реализация метода Рунге-Кутты для параллельного решения системы ОДУ»

textual
Листинг программы 
  1.    class Program
  2.     {
  3.         public delegate double Function(double x, double y);
  4.  
  5.         #region Ordinary Differential Equations - Methods
  6.  
  7.         public static double ODE_Euler(Function f, double x0, double y0, double h, double x)
  8.         {
  9.             double xnew, ynew, result = double.NaN;
  10.             if (x <= x0)
  11.                 result = y0;
  12.             else if (x > x0)
  13.             {
  14.                 do
  15.                 {
  16.                     if (h > x - x0) h = x - x0;
  17.                     ynew = y0 + f(x0, y0) * h;
  18.                     xnew = x0 + h;
  19.                     x0 = xnew;
  20.                     y0 = ynew;
  21.                 } while (x0 < x);
  22.                 result = ynew;
  23.             }
  24.             return result;
  25.         }
  26.  
  27.         public static double ODE_RungeKutta2(Function f, double x0, double y0, double h, double x)
  28.         {
  29.             double xnew, ynew, k1, k2, result = double.NaN;
  30.             if (x == x0)
  31.                 result = y0;
  32.             else if (x > x0)
  33.             {
  34.                 do
  35.                 {
  36.                     if (h > x - x0) h = x - x0;
  37.                     k1 = h * f(x0, y0);
  38.                     k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
  39.                     ynew = y0 + k2;
  40.                     xnew = x0 + h;
  41.                     x0 = xnew;
  42.                     y0 = ynew;
  43.                 } while (x0 < x);
  44.                 result = ynew;
  45.             }
  46.             return result;
  47.         }
  48.  
  49.         public static double ODE_RungeKutta4(Function f, double x0, double y0, double h, double x)
  50.         {
  51.             double xnew, ynew, k1, k2, k3, k4, result = double.NaN;
  52.             if (x == x0)
  53.                 result = y0;
  54.             else if (x > x0)
  55.             {
  56.                 do
  57.                 {
  58.                     if (h > x - x0) h = x - x0;
  59.                     k1 = h * f(x0, y0);
  60.                     k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
  61.                     k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2);
  62.                     k4 = h * f(x0 + h, y0 + k3);
  63.                     ynew = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
  64.                     xnew = x0 + h;
  65.                     x0 = xnew;
  66.                     y0 = ynew;                                        
  67.                 }while (x0 < x);
  68.                 result = ynew;
  69.             }
  70.             return result;
  71.         }
  72.  
  73.         public static double ODE_RungeKuttaFehlberg(Function f, double x0, double y0, double x, double h, double tolerance)
  74.         {
  75.             double xnew, ynew, hnew, k1, k2, k3, k4, k5, k6;
  76.             double hmin = 0.0001;
  77.             double hmax = 0.5;
  78.             if (h > hmax) h = hmax;
  79.             if (h < hmin) h = hmin;
  80.  
  81.             while (x0 < x)
  82.             {
  83.                 k1 = h * f(x0, y0);
  84.                 k2 = h * f(x0 + 0.25 * h, y0 + 0.25 * k1);
  85.                 k3 = h * f(x0 + 3 * h / 8, y0 + 3 * k1 / 32 + 9 * k2 / 32);
  86.                 k4 = h * f(x0 + 12 * h / 13, y0 + 1932 * k1 / 2197 - 7200 * k2 / 2197 + 7296 * k3 / 2197);
  87.                 k5 = h * f(x0 + h, y0 + 439 * k1 / 216 - 8 * k2 + 3680 * k3 / 513 - 845 * k4 / 4104);
  88.                 k6 = h * f(x0 + 0.5 * h, y0 - 8 * k1 / 27 + 2 * k2 - 3544 * k3 / 2565 + 1859 * k4 / 4104 - 11 * k5 / 40);
  89.                 double error = Math.Abs(k1 / 360 - 128 * k3 / 4275 - 2197 * k4 / 75240 + k5 / 50 + 2 * k6 / 55) / h;
  90.                 double s = Math.Pow(0.5 * tolerance / error, 0.25);
  91.                 if (error < tolerance)
  92.                 {
  93.                     ynew = y0 + 25 * k1 / 216 + 1408 * k3 / 2565 + 2197 * k4 / 4104 - 0.2 * k5;
  94.                     xnew = x0 + h;
  95.                     x0 = xnew;
  96.                     y0 = ynew; 
  97.                 }
  98.                 if (s < 0.1) s = 0.1;
  99.                 if (s > 4)   s = 4;
  100.                 hnew = h*s;
  101.                 h = hnew;
  102.                 if (h > hmax) h = hmax;
  103.                 if (h < hmin) h = hmin;
  104.                 if (h > x - x0) h = x - x0;
  105.             } return y0;
  106.         }
  107.  
  108.         #endregion
  109.  
  110.         #region Ordinary Differential Equations - Test functions
  111.  
  112.         static double f(double x, double y)
  113.         {
  114.             return y*Math.Cos(x);
  115.         }
  116.  
  117.         static double dx(double x, double y, double z)
  118.         { return 10.0 * (y - x); }
  119.  
  120.         static double dy(double x, double y, double z)
  121.         { return x * (28.0 - z) - y; }
  122.  
  123.         static double dz(double x, double y, double z)
  124.         { return x * y - 8.0 * z / 3.0; }
  125.  
  126.         static void TestEuler()
  127.         {
  128.             double h = 0.001;
  129.             double x0 = 0.0;
  130.             double y0 = 1.0;
  131.             Console.WriteLine("\n Results from Euler's method with h = {0}\n", h);
  132.             double result = y0;
  133.             for (int i = 0; i < 11; i++)
  134.             {
  135.                 double x = 0.1 * i;
  136.                 result = ODE_Euler(f, x0, result, h, x);
  137.                 double exact = Math.Exp(Math.Sin(x));
  138.                 if (i % 5 == 0)
  139.                 Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
  140.                 x0 = x;
  141.             }
  142.         }
  143.         
  144.         static void TestRungeKutta2()
  145.         {
  146.             double h = 0.001;
  147.             double x0 = 0.0;
  148.             double y0 = 1.0;
  149.             Console.WriteLine("\n Results from the 2nd-order Runge-Kutta method with h = {0}\n", h);
  150.             double result = y0;
  151.             for (int i = 0; i < 11; i++)
  152.             {
  153.                 double x = 0.1 * i;
  154.                 result = ODE_RungeKutta2(f, x0, result, h, x);
  155.                 double exact = Math.Exp(Math.Sin(x));
  156.                 if (i % 5 == 0)
  157.                 Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
  158.                 x0 = x;
  159.             }
  160.         }
  161.  
  162.         static void TestRungeKutta4()
  163.         {
  164.             double h = 0.001;
  165.             double x0 = 0.0;
  166.             double y0 = 1.0;
  167.             Console.WriteLine("\n Results from the 4th-order Runge-Kutta method with h = {0}\n", h);
  168.             double result = y0;
  169.             for (int i = 0; i < 11; i++)
  170.             {
  171.                 double x = 0.1 * i;
  172.                 result = ODE_RungeKutta4(f, x0, result, h, x);
  173.                 double exact = Math.Exp(Math.Sin(x));
  174.                 if (i % 5 == 0)
  175.                 Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
  176.                 x0 = x;
  177.             }
  178.         }
  179.  
  180.         static void TestRungeKuttaFehlberg()
  181.         {
  182.             double h = 0.2;
  183.             double x0 = 0.0;
  184.             double y0 = 1.0;
  185.             Console.WriteLine("\n Results from the fourth-order Runge-Kutta-Fehlberg method with h = {0}\n", h);
  186.             double result = y0;
  187.             for (int i = 0; i < 11; i++)
  188.             {
  189.                 double x = 0.1 * i;
  190.                 result = ODE_RungeKuttaFehlberg(f, x0, result, x, h, 1e-8);
  191.                 double exact = Math.Exp(Math.Sin(x));
  192.                 if (i%5==0)
  193.                 Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
  194.                 x0 = x;
  195.             }
  196.         }
  197.  
  198.         static void TestEulerCoupledODE()
  199.         {
  200.             double x = 0.0, y = 5.0, z = 10.0;  //initial conditions
  201.             double dt = 0.001;                  //step size
  202.  
  203.             Console.WriteLine("\n Results from a coupled ODE using Euler's method\n");
  204.  
  205.             for (int i = 0; i < 100; i++)
  206.             {
  207.                 double xnew = x + dx(x, y, z) * dt;
  208.                 double ynew = y + dy(x, y, z) * dt;
  209.                 double znew = z + dz(x, y, z) * dt;
  210.                 x = xnew;
  211.                 y = ynew;
  212.                 z = znew;
  213.                 Console.WriteLine(" x = {0:n1}, y = {1:e12}, z = {2:e12}", x, y, z);
  214.             }
  215.         }
  216.         #endregion
  217.         
  218.         static void Main(string[] args)
  219.         {
  220.             TestEuler();
  221.             TestRungeKutta2();
  222.             TestRungeKutta4();
  223.             TestRungeKuttaFehlberg();
  224.             TestEulerCoupledODE();
  225.             Console.ReadLine();
  226.         }
  227.     }

ИИ поможет Вам:

  • решить любую задачу по программированию
  • объяснить код
  • расставить комментарии в коде
  • и т.д
Попробуйте бесплатно

Оцени полезность:

8   голосов , оценка 4.375 из 5

Нужна аналогичная работа?

Оформи быстрый заказ и узнай стоимость

Бесплатно
Оформите заказ и авторы начнут откликаться уже через 10 минут
Похожие ответы