Решение системы линейных и дифференциальных уравнений - C#
Формулировка задачи:
Всем доброго времени суток.
Очень прошу помощи, ибо уже совсем отчаялся.
Необходимо решить систему из 12 уравнений - трёх дифференциальных и девяти алгебраических - на определённом временом отрезке.
Научным руководителем было предложено линеаризовать дифференциальные уравнения методом трапеций и затем решать полученную систему нелинейных алгебраических уравнений методом Ньютона, на каждом шаге решая её методом Гаусса.
У меня получилось записать всё это чудо в программу и даже получить адекватные числа. Одна проблема - я пишу замену для существующей досовской программы, и по идее данные должны с ней совпадать. Но они не совпадают, причём довольно сильно.
В связи с этим вопрос - возможно ли такой подход к решению?
И если да, то есть ли где-нибудь образцовая задача (не обязательно с таким количеством уравнений, можно хотя бы парочку ДУ и парочку АУ и правильные числа) для проверки правильности написанного кода?
Очень сильно прошу помощи. В гугле провёл долго, не нашёл.
Все переменные, не находящиеся слева в уравнениях - известные константы.
Решение задачи: «Решение системы линейных и дифференциальных уравнений»
textual
Листинг программы
class Program
{
public delegate double Function(double x, double y);
#region Ordinary Differential Equations - Methods
public static double ODE_Euler(Function f, double x0, double y0, double h, double x)
{
double xnew, ynew, result = double.NaN;
if (x <= x0)
result = y0;
else if (x > x0)
{
do
{
if (h > x - x0) h = x - x0;
ynew = y0 + f(x0, y0) * h;
xnew = x0 + h;
x0 = xnew;
y0 = ynew;
} while (x0 < x);
result = ynew;
}
return result;
}
public static double ODE_RungeKutta2(Function f, double x0, double y0, double h, double x)
{
double xnew, ynew, k1, k2, result = double.NaN;
if (x == x0)
result = y0;
else if (x > x0)
{
do
{
if (h > x - x0) h = x - x0;
k1 = h * f(x0, y0);
k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
ynew = y0 + k2;
xnew = x0 + h;
x0 = xnew;
y0 = ynew;
} while (x0 < x);
result = ynew;
}
return result;
}
public static double ODE_RungeKutta4(Function f, double x0, double y0, double h, double x)
{
double xnew, ynew, k1, k2, k3, k4, result = double.NaN;
if (x == x0)
result = y0;
else if (x > x0)
{
do
{
if (h > x - x0) h = x - x0;
k1 = h * f(x0, y0);
k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2);
k4 = h * f(x0 + h, y0 + k3);
ynew = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
xnew = x0 + h;
x0 = xnew;
y0 = ynew;
}while (x0 < x);
result = ynew;
}
return result;
}
public static double ODE_RungeKuttaFehlberg(Function f, double x0, double y0, double x, double h, double tolerance)
{
double xnew, ynew, hnew, k1, k2, k3, k4, k5, k6;
double hmin = 0.0001;
double hmax = 0.5;
if (h > hmax) h = hmax;
if (h < hmin) h = hmin;
while (x0 < x)
{
k1 = h * f(x0, y0);
k2 = h * f(x0 + 0.25 * h, y0 + 0.25 * k1);
k3 = h * f(x0 + 3 * h / 8, y0 + 3 * k1 / 32 + 9 * k2 / 32);
k4 = h * f(x0 + 12 * h / 13, y0 + 1932 * k1 / 2197 - 7200 * k2 / 2197 + 7296 * k3 / 2197);
k5 = h * f(x0 + h, y0 + 439 * k1 / 216 - 8 * k2 + 3680 * k3 / 513 - 845 * k4 / 4104);
k6 = h * f(x0 + 0.5 * h, y0 - 8 * k1 / 27 + 2 * k2 - 3544 * k3 / 2565 + 1859 * k4 / 4104 - 11 * k5 / 40);
double error = Math.Abs(k1 / 360 - 128 * k3 / 4275 - 2197 * k4 / 75240 + k5 / 50 + 2 * k6 / 55) / h;
double s = Math.Pow(0.5 * tolerance / error, 0.25);
if (error < tolerance)
{
ynew = y0 + 25 * k1 / 216 + 1408 * k3 / 2565 + 2197 * k4 / 4104 - 0.2 * k5;
xnew = x0 + h;
x0 = xnew;
y0 = ynew;
}
if (s < 0.1) s = 0.1;
if (s > 4) s = 4;
hnew = h*s;
h = hnew;
if (h > hmax) h = hmax;
if (h < hmin) h = hmin;
if (h > x - x0) h = x - x0;
} return y0;
}
#endregion
#region Ordinary Differential Equations - Test functions
static double f(double x, double y)
{
return y*Math.Cos(x);
}
static double dx(double x, double y, double z)
{ return 10.0 * (y - x); }
static double dy(double x, double y, double z)
{ return x * (28.0 - z) - y; }
static double dz(double x, double y, double z)
{ return x * y - 8.0 * z / 3.0; }
static void TestEuler()
{
double h = 0.001;
double x0 = 0.0;
double y0 = 1.0;
Console.WriteLine("\n Results from Euler's method with h = {0}\n", h);
double result = y0;
for (int i = 0; i < 11; i++)
{
double x = 0.1 * i;
result = ODE_Euler(f, x0, result, h, x);
double exact = Math.Exp(Math.Sin(x));
if (i % 5 == 0)
Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
x0 = x;
}
}
static void TestRungeKutta2()
{
double h = 0.001;
double x0 = 0.0;
double y0 = 1.0;
Console.WriteLine("\n Results from the 2nd-order Runge-Kutta method with h = {0}\n", h);
double result = y0;
for (int i = 0; i < 11; i++)
{
double x = 0.1 * i;
result = ODE_RungeKutta2(f, x0, result, h, x);
double exact = Math.Exp(Math.Sin(x));
if (i % 5 == 0)
Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
x0 = x;
}
}
static void TestRungeKutta4()
{
double h = 0.001;
double x0 = 0.0;
double y0 = 1.0;
Console.WriteLine("\n Results from the 4th-order Runge-Kutta method with h = {0}\n", h);
double result = y0;
for (int i = 0; i < 11; i++)
{
double x = 0.1 * i;
result = ODE_RungeKutta4(f, x0, result, h, x);
double exact = Math.Exp(Math.Sin(x));
if (i % 5 == 0)
Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
x0 = x;
}
}
static void TestRungeKuttaFehlberg()
{
double h = 0.2;
double x0 = 0.0;
double y0 = 1.0;
Console.WriteLine("\n Results from the fourth-order Runge-Kutta-Fehlberg method with h = {0}\n", h);
double result = y0;
for (int i = 0; i < 11; i++)
{
double x = 0.1 * i;
result = ODE_RungeKuttaFehlberg(f, x0, result, x, h, 1e-8);
double exact = Math.Exp(Math.Sin(x));
if (i%5==0)
Console.WriteLine(" x = {0:n1}, y = {1:e12}, exact = {2:e12}", x, result, exact);
x0 = x;
}
}
static void TestEulerCoupledODE()
{
double x = 0.0, y = 5.0, z = 10.0; //initial conditions
double dt = 0.001; //step size
Console.WriteLine("\n Results from a coupled ODE using Euler's method\n");
for (int i = 0; i < 100; i++)
{
double xnew = x + dx(x, y, z) * dt;
double ynew = y + dy(x, y, z) * dt;
double znew = z + dz(x, y, z) * dt;
x = xnew;
y = ynew;
z = znew;
Console.WriteLine(" x = {0:n1}, y = {1:e12}, z = {2:e12}", x, y, z);
}
}
#endregion
static void Main(string[] args)
{
TestEuler();
TestRungeKutta2();
TestRungeKutta4();
TestRungeKuttaFehlberg();
TestEulerCoupledODE();
Console.ReadLine();
}
ИИ поможет Вам:
- решить любую задачу по программированию
- объяснить код
- расставить комментарии в коде
- и т.д