Современные численные методы в объектно-ориентированном изложении на C#

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

Разбить на страницы
Показывать лекцию целиком

Цель лекции: Рассмотреть методы построения приближенных решений для линейных и нелинейных эволюционных уравнений.

В предыдущей лекции мы рассматривали системы обыкновенных дифференциальных уравнений конечной размерности. Однако многие процессы в нашем мире описываются бесконечными системами дифференциальных уравнений. Такие системы иногда называют распределенными системами или уравнениями в частных производных. Решением уравнений в частных производных является функция многих переменных. Простейшими уравнениями в частных производных второго порядка являются следующие уравнения$$\Delta u(x)=f(x)\quad\mbox{- уравнение Лапласа}$$ $$u_t(t,x)=\Delta u(t,x)+f(t,x)\quad\mbox{- уравнение теплопроводности}$$ $$u_{tt}(t,x)=\Delta u(t,x)+f(t,x)\quad\mbox{- волновое уравнение}$$ Уравнение Лапласа является примером эллиптического уравнения, уравнение теплопроводности является примером параболического уравнения, а волновое уравнение является примером гиперболического уравнения. С математической точки зрения дифференциальные уравнения в частных производных представляют собой весьма сложную тему, которая имеет принципиальные отличия от обыкновенных дифференциальных уравнений. С вычислительной точки зрения, часто, дифференциальные уравнения в частных производных могут быть аппроксимированы конечномерными уравнениями. Эллиптические уравнения могут быть аппроксимированы системами линейных алгебраических уравнений, а эволюционные уравнения аппроксимируются системами обыкновенных дифференциальных уравнений.

Мы будем рассматривать эволюционные уравнения относительно функций заданных на отрезке $$[0,T]$$ со значениями в банаховом пространстве, например, в пространстве $$C[a,b]$$. Численные методы, которые мы будем рассматривать могут быть одинаково эффективно применены как для линейных, так и для нелинейных уравнений. Рассмотрим математическую постановку начальной задачи. Пусть $$X$$ есть некоторое банахово пространство. Пусть в этом в этом пространстве задано множество $${\cal D}\subset X$$, на котором определен оператор $${\cal A}$$ (линейный или нелинейный).

Будем рассматривать уравнение$$u_t(t)={\cal A} u(t),\quad t\in[0,T]$$ с начальным условием$$u(0)=\varphi.$$ Искомой функцией в задаче 16.5-16.6 является функция $$u\in C^1([0,T];X)$$ такая, что $$u(t)\in{\cal D}$$ при всех $$t\in[0,T]$$ и, удовлетворяющая 16.5-16.6. Элемент $$\varphi\in{\cal D}$$ называется начальным условием, а задача 16.5-16.6 называется абстрактной задачей Коши.

Приведем два характерных примера таких задач. В качестве пространства $$X$$ мы возьмем пространство непрерывных функций $$C[0,2\pi]$$ в качестве множества $${\cal D}$$ возьмем множество непрерывно дифференцируемых функций на отрезке $$[0,2\pi]$$ и являющихся $$2\pi$$ -периодическими. Будем рассматривать два модельных уравнения. Первое - уравнение линейного переноса$$u_t(t,x)=cu_x(t,x),$$ где $$c\neq0$$ - постоянная, имеющая смысл "скорости" распространения волны. Второе уравнение - уравнение нелинейного переноса$$u_t(t,x)=u(t,x)u_x(t,x).$$

Опишем построение численной схемы. Мы будем использовать проекционный метод для приближенного решения абстрактной задачи Коши. Этот метод позволяет свести задачу к системе обыкновенных дифференциальных уравнений. Пусть для любого $$n>0$$ существует пара отображений$$P_N:X\to\Bbb{R}^n,\quad I_N:\Bbb{R}^n\to {\cal D}.$$ Как правило на эти отображения накладываются условия$$P_nI_nx=x,\quad\mbox{для любого\ }x\in\Bbb{R}^n$$ и$$\lim\limits_{n\to\infty}I_NP_Nx=x,\quad\mbox{для любого\ }x\in{\cal D},$$ где предел понимается в метрике пространства $$X$$. Задача 16.5-16.6 заменяется следующей задачей$$u^n_t(t)=P_n{\cal A} I_nu^n(t),$$ $$u^n_t(0)=P_n\\varphi.$$ Задача 16.9-16.10 представляет собой задачу Коши для системы обыкновенных дифференциальных уравнений $$n$$ -го порядка. Далее эта система решается стандартными численными методами, например, методом Рунге-Кутта, который мы рассматривали на прошлой лекции. После нахождения решения задачи 16.9-16.10 то есть функции $$u^n(t)$$, в качестве приближенным решением исходной задачи можно выбрать функцию $$I_nu^n(t)$$.

Хотя описанный проекционный метод является, как правило, легко реализуемым, но при его использовании необходимо иметь в виду вопросы, связанные с его сходимостью и устойчивостью. Дело в том, что даже для простейших уравнений при использовании проекционного метода следует согласовывать шаг по времени, то есть тот шаг, который используется в численном методе при решении задачи 16.9-16.10, с шагом, который имеет место при построении аппроксимации пространства $$X$$.

Для уравнений 16.7 и 16.8 мы будем использовать следующие операторы $$P_n$$ и $$I_n$$$$P_nf(x)=\left(% \begin{array}{c} f_1 \\ f_2 \\ \vdots \\ f_n \\ \end{array}% \right)$$ где $$f_i=f((i-1)\frac{2\pi}{n})$$, $$i=1,2,\dots,n$$. Для реализации оператора $$I_n$$ необходимо использовать подходящую интерполяцию. Мы будем рассматривать кусочно-линейную интерполяцию. Покажем, как таким образом можно определить операцию$$P_n\frac{d}{dx}I_nf=g=\left(% \begin{array}{c} g_1 \\ g_2 \\ \vdots \\ g_n \\ \end{array}% \right)$$ где$$g_i=\frac{f_{i+1}-f_i}{h},\quad i=1,2\dots,n-1,$$ $$g_n=\frac{f_{1}-f_n}{h},$$ где $$h=\frac{2\pi}{n}$$.

Реализуем этот метод для линейного уравнения переноса. Для этого мы создадим класс, являющийся наследником от класса $$TRungeKutta$$.

$$\begin{verbatim} class TCUEvol : TRungeKutta { double c; // скорость волны double h; public TCUEvol(int N, double c) : base(N) { this.c = c; h = 2.0 * Math.PI / N; } public override void F(double t, double[] Y, ref double[] FY) { int i; for (i = 0; i < N - 1; i++) { FY[i] = c * (Y[i + 1] - Y[i]) / h; } FY[N - 1] = c * (Y[0] - Y[N - 1]) / h; } } \end{verbatim}$$

Испытаем наш класс, учитывая, что для начальной функции вида $$\sin kx$$ это уравнение имеет точное решение в виде бегущей волны $$u(t,x)= sin(k(t+x))$$.

$$\begin{verbatim} int N = 1000; double h = 0.0001; double[] Y0 = new double[N]; TCUEvol CU = new TCUEvol(N, 1); double x; int i; for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Y0[i] = Math.Sin(5.0 * x); } CU.SetInit(0, Y0); double t = 0; while (CU.GetCurrent() < (1.0 + h / 2.0)) { t = CU.GetCurrent(); CU.NextStep(h); // рассчитать на следующем шаге } StreamWriter Fout = File.CreateText("cu.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}\t{2}\t{3}", x, CU.Y[i], Math.Sin(5.0 * (t + x)), Math.Sin(5.0 * (t + x)) - CU.Y[i]); } Fout.Close(); \end{verbatim}$$

Мы решали задачу:$$u_t(t,x)=u_x(t,x)$$ $$u(0,x)=\sin 5x.$$ Эта задача с $$2\pi$$ -периодическими по переменной $$x$$ условиями имеет единственное решение$$u(t,x)=\sin(5(t+x)),$$

в чем можно убедится непосредственно. На рисунке 17.1 мы приводим графики точного решения (сплошной линией) и приближенного (точечной линией) при $$t=1.0$$. Результат нашего вычислительного опыта не является удовлетворительным. Мы видим хорошее совпадение фазы решений, но амплитуда приближенного решения меньше точного. Наша схема оказалась весьма диссипативной, что делает ее непригодной. На рисунке 17.2 мы приведем погрешность приближенного решения.

(рис 17.2) Решение уравнения линейного переноса с помощью кусочно-линейной интерполяции(рис 17.1) Погрешность приближенного решения уравнения линейного переноса с помощью кусочно-линейной интерполяции

Покажем как можно существенно повысить точность нашего метода. Для этого мы воспользуемся так называемыми аналитико-числовыми методами. Суть этих методов состоит в том, что ряд операций, можно выполнить точно. Например в известной нам уже задаче дифференцирования функции. Пусть нам нужно найти производную функции $$f(x)=a\sin\alpha x+b\cos\beta x$$. При программировании нам выгодно задавать не поточечные значения этой функции, а всего четыре коэффициента - $$a$$, $$b$$, $$\alpha$$ и $$\beta$$. Обозначим эти коэффициенты четырьмя переменными:

$$\begin{verbatim} double sin_a; \\ a double sin_alpha; \\ alpha double cos_b; \\ b double cos_beta; \\ beta \end{verbatim}$$

Тогда производная функции $$f'(x)$$ имеет такой же вид и может быть также представлена такими же коэффициентами:

$$\begin{verbatim} double Dsin_a; \\ a double Dsin_alpha; \\ alpha double Dcos_b; \\ b double Dcos_beta; \\ beta \end{verbatim}$$

Вычислить коэффициенты производной можно следующим образом

$$\begin{verbatim} Dsin_a = -cos_b * cos_beta; Dsin_alpha = cos_beta; Dcos_b = sin_a * sin_alpha; Dcos_beta = sin_alpha; \end{verbatim}$$

При это, если не считать возможные ошибки округления при умножении, операция дифференцирования была выполнена безукоризненно.

С другой стороны гладкие (например, кусочно дифференцируемые и имеющие ограниченные производные) периодические функции могут быть разложены в ряд Фурье, а ряд Фурье представляет собой комбинацию рассмотренных выше функций. Поэтому, если задавать функцию коэффициентами Фурье, то такие операции как дифференцирование и интегрирование можно выполнять фантастически быстро. Чем более гладкая функция, тем быстрее коэффициенты Фурье сходятся к нулю. Следовательно, часто ряд Фурье удается аппроксимировать конечными тригонометрическими суммами.

Реализуем новый класс в котором мы будем применять численно-аналитические методы.

$$\begin{verbatim} class TCUAn : TRungeKutta { double c; double h; public TCUAn(int N, double c) : base(N) { this.c = c; h = 2.0 * Math.PI / N; } public double GetY(double x) { double res = 0; int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; res += Math.Sin(di * x) * Y[i - 1]; res += Math.Cos(di * x) * Y[N - i]; } res += Y[(N - 1) / 2]; return res; } \end{verbatim}$$ $$\begin{verbatim} public override void F(double t, double[] Y, ref double[] FY) { Diff(Y, ref FY); int i; for (i = 0; i < N; i++) { FY[i] = c * FY[i]; } } void Diff(double[] Y, ref double[] DY) { int i; double di; for (i = 1; i < (N-1) / 2; i++) { di = (double)i; DY[i - 1] = -di * Y[N - i]; DY[N - i] = di * Y[i - 1]; } DY[(N - 1) / 2] = 0; } } \end{verbatim}$$

В этом классе мы добавили новые методы - для дифференцирования функции, а также для суммирования тригонометрической суммы. Заметим, что мы рассматриваем систему дифференциальных уравнений относительно коэффициентов Фурье. Испытаем наш класс.

$$\begin{verbatim} int N = 257; double h = 0.001; double[] Y0 = new double[N]; TCUAn CUAn = new TCUAn(N, 1); for (i = 0; i < N; i++) { Y0[i] = 0; } Y0[4] = 1; CUAn.SetInit(0, Y0); t = 0; while (CUAn.GetCurrent() < (1.0 + h / 2.0)) { t = CUAn.GetCurrent(); CUAn.NextStep(h); } Fout = File.CreateText("cuan.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}\t{2}\t{3}", x, CUAn.GetY(x), Math.Sin(5.0 * (t + x)), Math.Sin(5.0 * (t + x)) - CUAn.GetY(x)); } Fout.Close(); \end{verbatim}$$

На рисунке 17.3 мы приводим погрешность приближенного решения.

(рис 17.3) Погрешность приближенного решения уравнения линейного переноса с помощью рядов Фурье

Сейчас мы рассмотрим нелинейное уравнения переноса$$u_t(t,x)=u(t,x)u_x(t,x)$$ также с периодическим краевыми условиями. Для этого уравнения характерен эффект обрушения волны. Действительно, в этом уравнении скорость переноса пропорциональна амплитуде. Поэтому "верхушка" волны движется быстрее основания, и в определенный момент производная решения будет стремится к бесконечности. Мы приведем приближенное решение уравнения нелинейного переноса. Для этого мы также будем использовать аналико-численные методы. Для этого нам нужно реализовать операцию произведения рядов Фурье, что вполне возможно. Приведем код этого класса.

$$\begin{verbatim} class TUU : TRungeKutta { double h; public TUU(int N) : base(N) { h = 2.0 * Math.PI / N; } public double GetY(double x) { double res = 0; int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; res += Math.Sin(di * x) * Y[i - 1]; res += Math.Cos(di * x) * Y[N - i]; } res += Y[(N - 1) / 2]; return res; } \end{verbatim}$$ $$\begin{verbatim} public override void F(double t, double[] Y, ref double[] FY) { double[] DU = new double[N]; Diff(Y, ref DU); Mult(Y, DU, ref FY); } int kSin(int k) { return k - 1; } int kCos(int k) { return N - k; } \end{verbatim}$$ $$\begin{verbatim} void Mult(double[] U, double[] DU, ref double[] UDU) { int i; for (i = 0; i < N; i++) { UDU[i] = 0; } int N2 = (N - 1)/2; int k, m; for (k = 1; k < (N - 1) / 2; k++) { for (m = 1; m < (N - 1) / 2; m++) { // sin kx * sin mx if (k != m) { UDU[kCos(Math.Abs(k - m))] += 0.5 * U[kSin(k)] * DU[kSin(m)]; } if (k + m < N2) { UDU[kCos(k + m)] += -0.5 * U[kSin(k)] * DU[kSin(m)]; } \end{verbatim}$$ $$\begin{verbatim} // sin kx * cos mx if (k + m < N2) { UDU[kSin(k + m)] += 0.5 * U[kSin(k)] * DU[kCos(m)]; } if (k > m) { UDU[kSin(k - m)] += 0.5 * U[kSin(k)] * DU[kCos(m)]; } if (k < m) { UDU[kSin(m - k)] += -0.5 * U[kSin(k)] * DU[kCos(m)]; } \end{verbatim}$$ $$\begin{verbatim} // cos kx * sin km if (k + m < N2) { UDU[kSin(k + m)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (m > k) { UDU[kSin(m - k)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (m < k) { UDU[kSin(k - m)] += -0.5 * U[kCos(k)] * DU[kSin(m)]; } // cos kx * cos mx if (k + m < N2) { UDU[kCos(k + m)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (k != m) { UDU[kCos(Math.Abs(k - m))] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } } } } \end{verbatim}$$ $$\begin{verbatim} void Diff(double[] Y, ref double[] DY) { int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; DY[i - 1] = -di * Y[N - i]; DY[N - i] = di * Y[i - 1]; } DY[(N - 1) / 2] = 0; } } \end{verbatim}$$

Теперь испытаем этот класс.

$$\begin{verbatim} int N = 257; double h = 0.001; double[] Y0 = new double[N]; TUU UU = new TUU(N); for (i = 0; i < N; i++) { Y0[i] = 0; } Y0[1] = 1; UU.SetInit(0, Y0); t = 0; while (UU.GetCurrent() < (0.3 + h / 2.0)) { t = UU.GetCurrent(); UU.NextStep(h); } Fout = File.CreateText("uu.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}", x, UU.GetY(x)); } Fout.Close(); \end{verbatim}$$

На графике 17.4 мы приведем график приближенного решения при $$t=0.3$$. К сожалению мы не имеем точного решения нелинейного уравнения, поэтому мы не приводим графика погрешности приближенного решения. Видно, что волна "пытается" обрушится.

(рис 17.4) Приближенное решение уравнения нелинейного переноса при t = 0.3

Ключевые термины

Распределенные системы - бесконечно мерные системы дифференциальных уравнений.

Начальное условие - функция, которой равно решение в начальный момент.

Абстрактная задача Коши - дифференциальное уравнение для функций со значениями в банаховом пространстве с заданным начальным условием.

Проекционный метод - метод аппроксимации эволюционных уравнений в частных производных системами обыкновенных дифференциальных уравнений.

Аналитико-числовые методы - численные методы которые содержат операции, выполняемые с использованием аналитического представления математических объектов.

Краткие итоги: Рассмотрены проекционные методы решения эволюционных уравнений в частных производных. Показано, что использование аналитико-числовых методов позволяет существенно улучшить точность расчетов.

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

Цель занятия

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

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

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

Мы уже отмечали, что задача Коши$$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}$$ Мы видим, что здесь имеет место пилообразное поведение траектории. Если мы будем решать эту задачу с меньшим шагом, то это приведет к тому, что амплитуда нашей "пилы" будет уменьшаться.

Лабораторная работа "Дифференциальные уравнения"

Цель занятия

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

Сценарий лабораторной работы

  • Рассмотреть задачу Коши$$y'(t)=y^2(t)+1$$ $$y(0)=0$$ которая не имеет решения на отрезке $$[0,T]$$, если $$T\ge\frac{\pi}{2}$$.
  • Написать программу на основе метода Рунге-Кутта.
  • Провести вычислительные опыты по приближенному решению этой задачи Коши на отрезках различной длинны.
  • Рассмотрите различные нелинейные уравнения в ограниченной области. С помощью вычислительных опытов найдите время выхода решения на границу области.
  • Для уравнений с известными решениями - сравните это время с результатами численных опытов.
    •

    Указания

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

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

    Страницы:

    Цель лекции: Рассмотреть методы построения приближенных решений для линейных и нелинейных эволюционных уравнений.

    В предыдущей лекции мы рассматривали системы обыкновенных дифференциальных уравнений конечной размерности. Однако многие процессы в нашем мире описываются бесконечными системами дифференциальных уравнений. Такие системы иногда называют распределенными системами или уравнениями в частных производных. Решением уравнений в частных производных является функция многих переменных. Простейшими уравнениями в частных производных второго порядка являются следующие уравнения$$\Delta u(x)=f(x)\quad\mbox{- уравнение Лапласа}$$ $$u_t(t,x)=\Delta u(t,x)+f(t,x)\quad\mbox{- уравнение теплопроводности}$$ $$u_{tt}(t,x)=\Delta u(t,x)+f(t,x)\quad\mbox{- волновое уравнение}$$ Уравнение Лапласа является примером эллиптического уравнения, уравнение теплопроводности является примером параболического уравнения, а волновое уравнение является примером гиперболического уравнения. С математической точки зрения дифференциальные уравнения в частных производных представляют собой весьма сложную тему, которая имеет принципиальные отличия от обыкновенных дифференциальных уравнений. С вычислительной точки зрения, часто, дифференциальные уравнения в частных производных могут быть аппроксимированы конечномерными уравнениями. Эллиптические уравнения могут быть аппроксимированы системами линейных алгебраических уравнений, а эволюционные уравнения аппроксимируются системами обыкновенных дифференциальных уравнений.

    Мы будем рассматривать эволюционные уравнения относительно функций заданных на отрезке $$[0,T]$$ со значениями в банаховом пространстве, например, в пространстве $$C[a,b]$$. Численные методы, которые мы будем рассматривать могут быть одинаково эффективно применены как для линейных, так и для нелинейных уравнений. Рассмотрим математическую постановку начальной задачи. Пусть $$X$$ есть некоторое банахово пространство. Пусть в этом в этом пространстве задано множество $${\cal D}\subset X$$, на котором определен оператор $${\cal A}$$ (линейный или нелинейный).

    Будем рассматривать уравнение$$u_t(t)={\cal A} u(t),\quad t\in[0,T]$$ с начальным условием$$u(0)=\varphi.$$ Искомой функцией в задаче 16.5-16.6 является функция $$u\in C^1([0,T];X)$$ такая, что $$u(t)\in{\cal D}$$ при всех $$t\in[0,T]$$ и, удовлетворяющая 16.5-16.6. Элемент $$\varphi\in{\cal D}$$ называется начальным условием, а задача 16.5-16.6 называется абстрактной задачей Коши.

    Приведем два характерных примера таких задач. В качестве пространства $$X$$ мы возьмем пространство непрерывных функций $$C[0,2\pi]$$ в качестве множества $${\cal D}$$ возьмем множество непрерывно дифференцируемых функций на отрезке $$[0,2\pi]$$ и являющихся $$2\pi$$ -периодическими. Будем рассматривать два модельных уравнения. Первое - уравнение линейного переноса$$u_t(t,x)=cu_x(t,x),$$ где $$c\neq0$$ - постоянная, имеющая смысл "скорости" распространения волны. Второе уравнение - уравнение нелинейного переноса$$u_t(t,x)=u(t,x)u_x(t,x).$$

    Опишем построение численной схемы. Мы будем использовать проекционный метод для приближенного решения абстрактной задачи Коши. Этот метод позволяет свести задачу к системе обыкновенных дифференциальных уравнений. Пусть для любого $$n>0$$ существует пара отображений$$P_N:X\to\Bbb{R}^n,\quad I_N:\Bbb{R}^n\to {\cal D}.$$ Как правило на эти отображения накладываются условия$$P_nI_nx=x,\quad\mbox{для любого\ }x\in\Bbb{R}^n$$ и$$\lim\limits_{n\to\infty}I_NP_Nx=x,\quad\mbox{для любого\ }x\in{\cal D},$$ где предел понимается в метрике пространства $$X$$. Задача 16.5-16.6 заменяется следующей задачей$$u^n_t(t)=P_n{\cal A} I_nu^n(t),$$ $$u^n_t(0)=P_n\\varphi.$$ Задача 16.9-16.10 представляет собой задачу Коши для системы обыкновенных дифференциальных уравнений $$n$$ -го порядка. Далее эта система решается стандартными численными методами, например, методом Рунге-Кутта, который мы рассматривали на прошлой лекции. После нахождения решения задачи 16.9-16.10 то есть функции $$u^n(t)$$, в качестве приближенным решением исходной задачи можно выбрать функцию $$I_nu^n(t)$$.

    Хотя описанный проекционный метод является, как правило, легко реализуемым, но при его использовании необходимо иметь в виду вопросы, связанные с его сходимостью и устойчивостью. Дело в том, что даже для простейших уравнений при использовании проекционного метода следует согласовывать шаг по времени, то есть тот шаг, который используется в численном методе при решении задачи 16.9-16.10, с шагом, который имеет место при построении аппроксимации пространства $$X$$.

    Для уравнений 16.7 и 16.8 мы будем использовать следующие операторы $$P_n$$ и $$I_n$$$$P_nf(x)=\left(% \begin{array}{c} f_1 \\ f_2 \\ \vdots \\ f_n \\ \end{array}% \right)$$ где $$f_i=f((i-1)\frac{2\pi}{n})$$, $$i=1,2,\dots,n$$. Для реализации оператора $$I_n$$ необходимо использовать подходящую интерполяцию. Мы будем рассматривать кусочно-линейную интерполяцию. Покажем, как таким образом можно определить операцию$$P_n\frac{d}{dx}I_nf=g=\left(% \begin{array}{c} g_1 \\ g_2 \\ \vdots \\ g_n \\ \end{array}% \right)$$ где$$g_i=\frac{f_{i+1}-f_i}{h},\quad i=1,2\dots,n-1,$$ $$g_n=\frac{f_{1}-f_n}{h},$$ где $$h=\frac{2\pi}{n}$$.

    Реализуем этот метод для линейного уравнения переноса. Для этого мы создадим класс, являющийся наследником от класса $$TRungeKutta$$.

    $$\begin{verbatim} class TCUEvol : TRungeKutta { double c; // скорость волны double h; public TCUEvol(int N, double c) : base(N) { this.c = c; h = 2.0 * Math.PI / N; } public override void F(double t, double[] Y, ref double[] FY) { int i; for (i = 0; i < N - 1; i++) { FY[i] = c * (Y[i + 1] - Y[i]) / h; } FY[N - 1] = c * (Y[0] - Y[N - 1]) / h; } } \end{verbatim}$$

    Испытаем наш класс, учитывая, что для начальной функции вида $$\sin kx$$ это уравнение имеет точное решение в виде бегущей волны $$u(t,x)= sin(k(t+x))$$.

    $$\begin{verbatim} int N = 1000; double h = 0.0001; double[] Y0 = new double[N]; TCUEvol CU = new TCUEvol(N, 1); double x; int i; for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Y0[i] = Math.Sin(5.0 * x); } CU.SetInit(0, Y0); double t = 0; while (CU.GetCurrent() < (1.0 + h / 2.0)) { t = CU.GetCurrent(); CU.NextStep(h); // рассчитать на следующем шаге } StreamWriter Fout = File.CreateText("cu.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}\t{2}\t{3}", x, CU.Y[i], Math.Sin(5.0 * (t + x)), Math.Sin(5.0 * (t + x)) - CU.Y[i]); } Fout.Close(); \end{verbatim}$$

    Мы решали задачу:$$u_t(t,x)=u_x(t,x)$$ $$u(0,x)=\sin 5x.$$ Эта задача с $$2\pi$$ -периодическими по переменной $$x$$ условиями имеет единственное решение$$u(t,x)=\sin(5(t+x)),$$

    в чем можно убедится непосредственно. На рисунке 17.1 мы приводим графики точного решения (сплошной линией) и приближенного (точечной линией) при $$t=1.0$$. Результат нашего вычислительного опыта не является удовлетворительным. Мы видим хорошее совпадение фазы решений, но амплитуда приближенного решения меньше точного. Наша схема оказалась весьма диссипативной, что делает ее непригодной. На рисунке 17.2 мы приведем погрешность приближенного решения.

    (рис 17.2) Решение уравнения линейного переноса с помощью кусочно-линейной интерполяции(рис 17.1) Погрешность приближенного решения уравнения линейного переноса с помощью кусочно-линейной интерполяции

    Покажем как можно существенно повысить точность нашего метода. Для этого мы воспользуемся так называемыми аналитико-числовыми методами. Суть этих методов состоит в том, что ряд операций, можно выполнить точно. Например в известной нам уже задаче дифференцирования функции. Пусть нам нужно найти производную функции $$f(x)=a\sin\alpha x+b\cos\beta x$$. При программировании нам выгодно задавать не поточечные значения этой функции, а всего четыре коэффициента - $$a$$, $$b$$, $$\alpha$$ и $$\beta$$. Обозначим эти коэффициенты четырьмя переменными:

    $$\begin{verbatim} double sin_a; \\ a double sin_alpha; \\ alpha double cos_b; \\ b double cos_beta; \\ beta \end{verbatim}$$

    Тогда производная функции $$f'(x)$$ имеет такой же вид и может быть также представлена такими же коэффициентами:

    $$\begin{verbatim} double Dsin_a; \\ a double Dsin_alpha; \\ alpha double Dcos_b; \\ b double Dcos_beta; \\ beta \end{verbatim}$$

    Вычислить коэффициенты производной можно следующим образом

    $$\begin{verbatim} Dsin_a = -cos_b * cos_beta; Dsin_alpha = cos_beta; Dcos_b = sin_a * sin_alpha; Dcos_beta = sin_alpha; \end{verbatim}$$

    При это, если не считать возможные ошибки округления при умножении, операция дифференцирования была выполнена безукоризненно.

    С другой стороны гладкие (например, кусочно дифференцируемые и имеющие ограниченные производные) периодические функции могут быть разложены в ряд Фурье, а ряд Фурье представляет собой комбинацию рассмотренных выше функций. Поэтому, если задавать функцию коэффициентами Фурье, то такие операции как дифференцирование и интегрирование можно выполнять фантастически быстро. Чем более гладкая функция, тем быстрее коэффициенты Фурье сходятся к нулю. Следовательно, часто ряд Фурье удается аппроксимировать конечными тригонометрическими суммами.

    Реализуем новый класс в котором мы будем применять численно-аналитические методы.

    $$\begin{verbatim} class TCUAn : TRungeKutta { double c; double h; public TCUAn(int N, double c) : base(N) { this.c = c; h = 2.0 * Math.PI / N; } public double GetY(double x) { double res = 0; int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; res += Math.Sin(di * x) * Y[i - 1]; res += Math.Cos(di * x) * Y[N - i]; } res += Y[(N - 1) / 2]; return res; } \end{verbatim}$$ $$\begin{verbatim} public override void F(double t, double[] Y, ref double[] FY) { Diff(Y, ref FY); int i; for (i = 0; i < N; i++) { FY[i] = c * FY[i]; } } void Diff(double[] Y, ref double[] DY) { int i; double di; for (i = 1; i < (N-1) / 2; i++) { di = (double)i; DY[i - 1] = -di * Y[N - i]; DY[N - i] = di * Y[i - 1]; } DY[(N - 1) / 2] = 0; } } \end{verbatim}$$

    В этом классе мы добавили новые методы - для дифференцирования функции, а также для суммирования тригонометрической суммы. Заметим, что мы рассматриваем систему дифференциальных уравнений относительно коэффициентов Фурье. Испытаем наш класс.

    $$\begin{verbatim} int N = 257; double h = 0.001; double[] Y0 = new double[N]; TCUAn CUAn = new TCUAn(N, 1); for (i = 0; i < N; i++) { Y0[i] = 0; } Y0[4] = 1; CUAn.SetInit(0, Y0); t = 0; while (CUAn.GetCurrent() < (1.0 + h / 2.0)) { t = CUAn.GetCurrent(); CUAn.NextStep(h); } Fout = File.CreateText("cuan.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}\t{2}\t{3}", x, CUAn.GetY(x), Math.Sin(5.0 * (t + x)), Math.Sin(5.0 * (t + x)) - CUAn.GetY(x)); } Fout.Close(); \end{verbatim}$$

    На рисунке 17.3 мы приводим погрешность приближенного решения.

    (рис 17.3) Погрешность приближенного решения уравнения линейного переноса с помощью рядов Фурье

    Сейчас мы рассмотрим нелинейное уравнения переноса$$u_t(t,x)=u(t,x)u_x(t,x)$$ также с периодическим краевыми условиями. Для этого уравнения характерен эффект обрушения волны. Действительно, в этом уравнении скорость переноса пропорциональна амплитуде. Поэтому "верхушка" волны движется быстрее основания, и в определенный момент производная решения будет стремится к бесконечности. Мы приведем приближенное решение уравнения нелинейного переноса. Для этого мы также будем использовать аналико-численные методы. Для этого нам нужно реализовать операцию произведения рядов Фурье, что вполне возможно. Приведем код этого класса.

    $$\begin{verbatim} class TUU : TRungeKutta { double h; public TUU(int N) : base(N) { h = 2.0 * Math.PI / N; } public double GetY(double x) { double res = 0; int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; res += Math.Sin(di * x) * Y[i - 1]; res += Math.Cos(di * x) * Y[N - i]; } res += Y[(N - 1) / 2]; return res; } \end{verbatim}$$ $$\begin{verbatim} public override void F(double t, double[] Y, ref double[] FY) { double[] DU = new double[N]; Diff(Y, ref DU); Mult(Y, DU, ref FY); } int kSin(int k) { return k - 1; } int kCos(int k) { return N - k; } \end{verbatim}$$ $$\begin{verbatim} void Mult(double[] U, double[] DU, ref double[] UDU) { int i; for (i = 0; i < N; i++) { UDU[i] = 0; } int N2 = (N - 1)/2; int k, m; for (k = 1; k < (N - 1) / 2; k++) { for (m = 1; m < (N - 1) / 2; m++) { // sin kx * sin mx if (k != m) { UDU[kCos(Math.Abs(k - m))] += 0.5 * U[kSin(k)] * DU[kSin(m)]; } if (k + m < N2) { UDU[kCos(k + m)] += -0.5 * U[kSin(k)] * DU[kSin(m)]; } \end{verbatim}$$ $$\begin{verbatim} // sin kx * cos mx if (k + m < N2) { UDU[kSin(k + m)] += 0.5 * U[kSin(k)] * DU[kCos(m)]; } if (k > m) { UDU[kSin(k - m)] += 0.5 * U[kSin(k)] * DU[kCos(m)]; } if (k < m) { UDU[kSin(m - k)] += -0.5 * U[kSin(k)] * DU[kCos(m)]; } \end{verbatim}$$ $$\begin{verbatim} // cos kx * sin km if (k + m < N2) { UDU[kSin(k + m)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (m > k) { UDU[kSin(m - k)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (m < k) { UDU[kSin(k - m)] += -0.5 * U[kCos(k)] * DU[kSin(m)]; } // cos kx * cos mx if (k + m < N2) { UDU[kCos(k + m)] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } if (k != m) { UDU[kCos(Math.Abs(k - m))] += 0.5 * U[kCos(k)] * DU[kSin(m)]; } } } } \end{verbatim}$$ $$\begin{verbatim} void Diff(double[] Y, ref double[] DY) { int i; double di; for (i = 1; i < (N - 1) / 2; i++) { di = (double)i; DY[i - 1] = -di * Y[N - i]; DY[N - i] = di * Y[i - 1]; } DY[(N - 1) / 2] = 0; } } \end{verbatim}$$

    Теперь испытаем этот класс.

    $$\begin{verbatim} int N = 257; double h = 0.001; double[] Y0 = new double[N]; TUU UU = new TUU(N); for (i = 0; i < N; i++) { Y0[i] = 0; } Y0[1] = 1; UU.SetInit(0, Y0); t = 0; while (UU.GetCurrent() < (0.3 + h / 2.0)) { t = UU.GetCurrent(); UU.NextStep(h); } Fout = File.CreateText("uu.txt"); for (i = 0; i < N; i++) { x = (double)i * (2.0 * Math.PI) / (double)N; Fout.WriteLine("{0}\t{1}", x, UU.GetY(x)); } Fout.Close(); \end{verbatim}$$

    На графике 17.4 мы приведем график приближенного решения при $$t=0.3$$. К сожалению мы не имеем точного решения нелинейного уравнения, поэтому мы не приводим графика погрешности приближенного решения. Видно, что волна "пытается" обрушится.

    (рис 17.4) Приближенное решение уравнения нелинейного переноса при t = 0.3

    Ключевые термины

    Распределенные системы - бесконечно мерные системы дифференциальных уравнений.

    Начальное условие - функция, которой равно решение в начальный момент.

    Абстрактная задача Коши - дифференциальное уравнение для функций со значениями в банаховом пространстве с заданным начальным условием.

    Проекционный метод - метод аппроксимации эволюционных уравнений в частных производных системами обыкновенных дифференциальных уравнений.

    Аналитико-числовые методы - численные методы которые содержат операции, выполняемые с использованием аналитического представления математических объектов.

    Краткие итоги: Рассмотрены проекционные методы решения эволюционных уравнений в частных производных. Показано, что использование аналитико-числовых методов позволяет существенно улучшить точность расчетов.

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

    Цель занятия

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

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

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

    Мы уже отмечали, что задача Коши$$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}$$ Мы видим, что здесь имеет место пилообразное поведение траектории. Если мы будем решать эту задачу с меньшим шагом, то это приведет к тому, что амплитуда нашей "пилы" будет уменьшаться.

    Лабораторная работа "Дифференциальные уравнения"

    Цель занятия

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

    Сценарий лабораторной работы

  • Рассмотреть задачу Коши$$y'(t)=y^2(t)+1$$ $$y(0)=0$$ которая не имеет решения на отрезке $$[0,T]$$, если $$T\ge\frac{\pi}{2}$$.
  • Написать программу на основе метода Рунге-Кутта.
  • Провести вычислительные опыты по приближенному решению этой задачи Коши на отрезках различной длинны.
  • Рассмотрите различные нелинейные уравнения в ограниченной области. С помощью вычислительных опытов найдите время выхода решения на границу области.
  • Для уравнений с известными решениями - сравните это время с результатами численных опытов.
    •

    Указания

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

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

    Предыдущее занятие Следующее занятие
    Вернуться к учебному плану
    ^