Опубликован: 18.05.2011 | Доступ: свободный | Студентов: 871 / 65 | Оценка: 4.40 / 4.20 | Длительность: 12:30:00
Лекция 18:

Эволюционные уравнения в частных производных

< Лекция 17 || Лекция 18: 12345 || Лекция 19 >

Практическое занятие "Дифференциальные уравнения"

Цель занятия

Выяснить как ведут себя численные методы решения обыкновенных дифференциальных уравнений в случаях, когда отсутствует решение задачи.

Практическая задача

В соответствующей лекции мы рассматривали метод Рунге-Кутта для построения приближенных решений задачи Коши для системы обыкновенных дифференциальных уравнений. Сейчас мы рассмотрим некоторые примеры численных опытов.

Мы уже отмечали, что задача Коши

y'(t)=y^2(t)+1
y(0)=0
не имеет решения на отрезке [0,T], где T\ge\frac{\pi}{2}. Посмотрим как себя поведет численный метод Эйлера при попытке решить эту задачу на отрезке [0,2]. Выполним следующий код, используя написанный ранее класс TEuler.

\begin{verbatim}
double h = 0.1;

double[] Y0 = { 0 };

TTest1Euler Test1Euler = new TTest1Euler();
Test1Euler.SetInit(0, Y0);

while (Test1Euler.GetCurrent() < (2.0 + h / 2.0))
{
    Console.WriteLine("{0}\t{1}", Test1Euler.GetCurrent(),
    Test1Euler.Y[0]);

    Test1Euler.NextStep(h);
}
\end{verbatim}

Здесь используется следующий класс TTest1Euler.

\begin{verbatim}
class TTest1Euler : TEuler
{
    public TTest1Euler() : base(1) { }
    public override void F(double t, double[] Y, ref double[] FY)
    {
        FY[0] = Y[0] * Y[0] + 1.0;
    }
}
\end{verbatim}

Сначала мы используем нарочито грубый шаг по времени. Поэтому вот результат:

\begin{verbatim}
0       0

0.1     0.1

0.2     0.201

0.3     0.3050401

0.4     0.414345046260801

0.5     0.531513227996888

0.6     0.659763859150455

0.7     0.803292694134565

0.8     0.967820609379562

0.9     1.16148828257354

1       1.39639378562911

1.1     1.69138534608347

1.2     2.07746378497806

1.3     2.60904936276759

1.4     3.38976322050339

1.5     4.63881268961114

1.6     6.89067100654088

1.7     11.7388056985792

1.8     25.6187616214787

1.9     91.3508563232936

2       925.948751423197
\end{verbatim}

Мы видим, что наш метод в шагом h=10^{-1} практически "не чувствует", что решения уже не существует. Теперь уменьшим шаг на порядок - возьмем h=10^{-2}. В результате получим следующее.

\begin{verbatim}
1.49    9.66252838596756

1.5     10.6061729340638

1.51    11.7410819771365

1.52    13.1296120370749

1.53    14.863479159516

1.54    17.0827092867696

1.55    20.0108988525325

1.56    24.0252595813953

1.57    29.8073905609296

1.58    38.7021958814475

1.59    53.6907955419068

1.6     82.5278108011353

1.61    150.646206357415
\end{verbatim}
\begin{verbatim}
1.62    377.599001256224

1.63    1803.4190587532

1.64    34326.6320734961

1.65    11817503.3371652

1.66    1396545668742.45

1.67    1.95033980502295E+22

1.68    3.80382535505695E+42

1.69    1.44690873317741E+83

1.7     2.09354488214506E+164

1.71    бесконечность

1.72    бесконечность
\end{verbatim}

Мы видим, что C#, точнее .NET Framework, довольно корректно использовали значение "бесконечность". В этом случае уже можно говорить, что наш метод "заметил" отсутствие решения.

Теперь мы рассмотрим поведение нашего численного метода в случае, когда имеет место разрыв правой части по фазовым переменным. В этом случае поведение решений может быть по разному. Возможны случаи, когда уравнение имеет разрывные траектории, а возможны и случаи, когда уравнение не имеет решений. В последнем случае используют аппроксимацию дифференциального уравнения дифференциальным включением. Однако мы рассмотрим применение нашего метода Эйлера "в лоб" для задачи Коши

y'(t)=1-2\sign y(t),
y(0)=0.
Можно показать, что эта задача не имеет решения. Однако мы можем применить к ней метод Эйлера. Приведем листинг нашего класса, который решает эту задачу.

\begin{verbatim}
class TTest2Euler : TEuler
{
    public TTest2Euler() : base(1) { }
    public override void F(double t, double[] Y,
    ref double[] FY)
    {
        FY[0] = 1.0 - 2 * sign(Y[0]);
    }
    double sign(double x)
    {
        if (x >= 0)
        {
            return 1;
        }
        else
        {
            return -1;
        }
    }
}
\end{verbatim}

Запустим наш класс с шагом h=10^{-1} на отрезке [0,1]. Вот результат.

\begin{verbatim}
0       0

0.1     -0.1

0.2     0.2

0.3     0.1

0.4     2.77555756156289E-17

0.5     -0.1

0.6     0.2

0.7     0.1

0.8     5.55111512312578E-17

0.9     -0.1

1       0.2
\end{verbatim}
Мы видим, что здесь имеет место пилообразное поведение траектории. Если мы будем решать эту задачу с меньшим шагом, то это приведет к тому, что амплитуда нашей "пилы" будет уменьшаться.

< Лекция 17 || Лекция 18: 12345 || Лекция 19 >
Олег Корсак
Олег Корсак
Латвия, Рига
Дмитрий Кифель
Дмитрий Кифель
Казахстан, Темиртау