Нужна помощь с методом конечных разностей

[kometa_triatlon]

Здравствуйте.
Помогите найти ошибку в программе, которая решает линейную краевую задачу методом конечных разностей.
Даны следующие условия:
$ \[ \begin{gathered} ay'' - y = cos^2 \pi x + 2a\pi ^2 cos2\pi x, \hfill \\ a = \frac{1} {{400}}, y(0) = 0,y(1) = 0 \hfill \\ \end{gathered} \] $

Отсюда получаем:
$\[ \begin{gathered} p(x) = 0; \hfill \\ q(x) = 400; \hfill \\ f(x) = 400\cos ^2 \pi x + 2\pi ^2 \cos 2\pi x; \hfill \\ \end{gathered} \] $

Код программы на С++ (решение СЛАУ я опустил):

#define PI 3.1415926535897

double p(double x)
{
        return 0;
}

double q(double x)
{
        return 400;
}

double f(double x)
{
        return 400 * pow(cos(PI*x), 2) + 2 * PI * PI *cos(2*PI*x);
}

void main()
{
        double a = 400;
        long n = 16;

        double x0 = 0;
        double xn = 1;
        double y0 = 0;
        double yn = 0;
        double alpha0 = 1;
        double beta0 = 1;
        double h = fabs(xn - x0) / (double)(n - 1);

        double* x = new double[n];
        for(int i = 0; i < n; i++)
            x[i] = x0 + i*h;

        double** A = new double*[n];  // Трехдиагональная матрица
        double* B = new double[n];     // Вектор правой части

        for(int i = 0; i < n ; i++)
        {
            A[i] = new double[n];
            for(int j = 0; j < n; j++)
                A[i][j] = 0;
        }

        A[0][0] = alpha0;
        B[0] =  y0;
        for(int i = 0; i < n-2; i++)
        {

            A[i + 1][i] = 1 - h*p(x[i]) + h*h*q(x[i]);
            A[i + 1][i + 1] = h*p(x[i]) - 2;
            A[i + 1][i + 2] = 1;
            B[i + 1] = f(x[i]) * h * h;
        }

        A[n-1][n-1] = beta0;
        B[n-1] = yn;
        
        /// ....
        // Решаем СЛАУ методом Гаусса
}

Вроде все работает, но результаты - фигня полная
Вот пример:

N = 4
Матрица А:

   1.000   0.000   0.000   0.000
  45.444  -2.000   1.000   0.000
   0.000  45.444  -2.000   1.000
   0.000   0.000   0.000   1.000

Вектор свободных коэффициентов:

   0.000
  46.638
  10.014
   0.000

Решение:

   0.000
   2.492
  51.622
   0.000

Вроде все хорошо, но при увеличении количества интервалов искомая функция ведет себя очень странно:

N = 6
Матрица A:

   1.000   0.000   0.000   0.000   0.000   0.000
  17.000  -2.000   1.000   0.000   0.000   0.000
   0.000  17.000  -2.000   1.000   0.000   0.000
   0.000   0.000  17.000  -2.000   1.000   0.000
   0.000   0.000   0.000  17.000  -2.000   1.000
   0.000   0.000   0.000   0.000   0.000   1.000

Вектор свободных коэффициентов:

  
  0.000
  16.790
  10.716
   0.889
   0.889
   0.000

Решение:

   0.000
  11.327
  39.443
-102.954
-875.557
   0.000

При N = 10 решение "улетает" уже в положительную сторону:

   0.000
  -0.300
   3.597
  12.473
   9.638
 -41.766
-131.500
 -54.367
 548.987
1371.128
   0.000

При дальнейшем увеличении N вроде как приходит точность:

0.000
  -1.756
  -3.250
  -4.045
  -3.771
  -2.238
   0.473
   3.966
   7.546
  10.323
  11.380
  10.003
   5.906
  -0.588
  -8.476
 -16.155
 -21.670
 -23.118
 -19.135
  -9.370
   5.173
  22.048
  37.618
  47.665
  48.301
  37.025
  13.686
 -18.880
 -54.824
 -85.985
-103.358
 -99.130
 -68.939
 -13.818
  58.705
 134.870
 196.564
 224.763
 204.057
 127.408
   0.000

Вобщем, я что-то сделал не так, но не пойму что
Буду благодарен за любую помощь, может у кого-то найдутся исходники?

P.S. Такую СЛАУ вообще-то надо методом прогонки решать, я знаю. Пробовал, тоже фигня получается, мне кажется у меня ошибки в заполнении матриц, поэтому этот момент я пока упускаю.

[Фантом]

kometa_triatlon в сообщении #273965 писал(а):

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

Ошибка не в программе (точнее, возможно, где-то в коде какие-то мелочи и есть, но это непринципиально). Ошибка в выборе метода решения задачи.

Судя по условию, задача учебная, причем Вам специально подложили грабли для того, чтобы Вы набили на них шишку. Эта затея удалась, с чем я Вашего преподавателя и поздравляю.

Дело вот в чем. Давайте посмотрим на характеристическое уравнение для соответствующего однородного дифференциального уравнения. Оно имеет вид $a\lambda^2 - 1 = 0$. Решения, соответственно, $\lambda_{1,2}=\pm 1/\sqrt{a}$. При имеющемся значении параметра это $\lambda_{1,2}=\pm 20$.

Это означает, что любое решение уравнения имеет вид $y=C_1\,e^{20\,x}+C_2\,e^{-20\,x}+R(x)$, причем функция $R(x)$, вообще говоря, является комбинацией чего-то тригонометрического (и, как следствие, не слишком сильно меняется с изменением $x$).

А теперь представьте себе, как будет вести себя любое численное решение Вашей краевой задачи из-за банальных вычислительных погрешностей. Для того, чтобы это решение, стартовав в $y(0)=0$, пришло в $y(1)=0$ (или наоборот, это непринципиально), коэффициенты $C_{1,2}$ должны быть найдены крайне точно, иначе экспоненты с большим по модулю показателем обеспечат вместо решения нечто, напоминающее генератор случайных чисел. Соответственно, выбранный Вами метод требует либо вычислений с огромной точностью (а не жалкие 15 знаков типа double), либо использования очень большого количества интервалов.

Для того, чтобы всех этих неприятностей не случилось, существует полезная вещь - метод дифференциальной прогонки. Вот о нем и почитайте (посмотрите в интернете).

[kometa_triatlon]

Вау. Спасибо огромное....
Я просто в шоке, буду искать.
Теперь я хоть понял что у меня происходит с программой. Я вообще-то не математик, а программист, таких тонкостей не знал.
Причем ведь не мехмат какой-нибудь, а второй курс факультета менеджмента и маркетинга. Наши преподы злые

-- Вт дек 22, 2009 01:41:03 --

Кстати, в условии задания сказано решить задачу методом конечных разностей второго порядка.
Тут у меня был первый порядок, но я так понимаю второй меня не спасет?
Видимо остается только мириться с такими погрешностями...

[Фантом]

kometa_triatlon в сообщении #273980 писал(а):

Я вообще-то не математик, а программист, таких тонкостей не знал.
Причем ведь не мехмат какой-нибудь, а второй курс факультета менеджмента и маркетинга.

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

kometa_triatlon в сообщении #273980 писал(а):

Кстати, в условии задания сказано решить задачу методом конечных разностей второго порядка.
Тут у меня был первый порядок, но я так понимаю второй меня не спасет?

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

[kometa_triatlon]

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

[ewert]

kometa_triatlon в сообщении #273965 писал(а):

N = 6
Матрица A:

   1.000   0.000   0.000   0.000   0.000   0.000
  17.000  -2.000   1.000   0.000   0.000   0.000
   0.000  17.000  -2.000   1.000   0.000   0.000
   0.000   0.000  17.000  -2.000   1.000   0.000
   0.000   0.000   0.000  17.000  -2.000   1.000
   0.000   0.000   0.000   0.000   0.000   1.000

Это какая-то загадка, а не матрица. Она ж должна быть симметричной! Причём с жутко доминирующей главной диагональю. Соответственно (при таком-то малом количестве шагов) -- жутко хорошо обусловленной, и как ни решай систему -- одинаково хороший результат получится.

Фантом в сообщении #273975 писал(а):

Это означает, что любое решение уравнения имеет вид $y=C_1\,e^{20\,x}+C_2\,e^{-20\,x}+R(x)$, причем функция $R(x)$, вообще говоря, является комбинацией чего-то тригонометрического (и, как следствие, не слишком сильно меняется с изменением $x$).

А теперь представьте себе, как будет вести себя любое численное решение Вашей краевой задачи из-за банальных вычислительных погрешностей.

Очень хорошо будет себя вести. Наличие таких крутых экспонент как раз и означает, что мы о-очень далеки от спектра. Задача-то не пристрелкой решается, а системой.

[Фантом]

ewert в сообщении #274004 писал(а):

Очень хорошо будет себя вести. Наличие таких крутых экспонент как раз и означает, что мы о-очень далеки от спектра. Задача-то не пристрелкой решается, а системой.

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

[kometa_triatlon]

Цитата:

Это какая-то загадка, а не матрица.

Вполне возможно
Я использовал правую аппроксимацию производной:
$\[\begin{gathered} y_i^' \approx \frac{{y_{i + 1} - y_i }} {h}, \hfill \\ y{}_i^{''} \approx \frac{{y_{i + 2}^{} - 2y_{i + 1} + y_i }} {{h^2 }} \hfill \\ \end{gathered} \]$

Система уравнений имеет вид:
$\[ (1 - hp_i + h^2 q_i )y_i + ( - 2 + hp_i )y_{i + 1} + y_{i + 2} = f_i \cdot h^2 \] $

Таким образом заполняются диагонали на строках со 2-й по n-1.
И еще краевые условия дают по одному уравнению (первая и последняя строка матрицы):
$\[y_0 = A\]\[y_n = B\]$

Функции p и q от х не зависят, поэтому в каждой строчке получаем те же числа.
Вроде все правильно, или я где-то ошибся?

Вчера эксперементировал с программой, получил следующие результаты:

• Если количество интервалов малое, решение действительно скачет как генератор рандома
• При N больше тысячи начинает сходиться к точному решению
• Параметр а у меня вводится пользователем, поигрался с ним, действительно: если а небольшое решение ведет себя "хорошо": с увеличением количества интервалов график не скачет, а сглаживается, уточняя найденную функцию

-- Вт дек 22, 2009 11:58:31 --

Бонус-трек: визуализация
При а = 400:






-- Вт дек 22, 2009 12:02:37 --

а = 8:

[ewert]

kometa_triatlon в сообщении #274049 писал(а):

Я использовал правую аппроксимацию производной:
$\[\begin{gathered} y_i^' \approx \frac{{y_{i + 1} - y_i }} {h}, \hfill \\ y{}_i^{''} \approx \frac{{y_{i + 2}^{} - 2y_{i + 1} + y_i }} {{h^2 }} \hfill \\ \end{gathered} \]$

Хрен знает, чему, как и зачем вас учили.

Используйте самую тупую, стандартную симметричную аппроксимацию для второй производной: $y''_i\approx{y_{i-1}-2y_{i}+y_{i+1}\over h^2}$. И получите прекрасную, симметричную (ну с точностью до граничных условий) и со всех точек зрения восхитительную матрицу.

А так -- Вы, судя по всему, попросту не ту диагональ и не к той прибавили. Хотя даже и в этом случае: цифирка 17 (и это при $n=6$-то) вызывает глубокие подозрения.

[kometa_triatlon]

Цитата:

Хрен знает, чему, как и зачем вас учили.

А меня не учили
Это мне самому появилась надобность разобраться. Нашел реферат, где эти схемы подробно описаны, реализовал. Теперь вот с вашей помощью разбираюсь.
Кстати, аппроксимация которую вы привели имеет порядок точности $O(h^2)$, а моя - $O(h)$, так что она еще не самая тупая))
Ок, спасибо, у меня и самого была мысля попробовать другой вариант аппроксимации.

[ewert]

Фантом в сообщении #274040 писал(а):

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

Не смогут. Поскольку обратная матрица существенно ограничена -- задача устойчива. И, следовательно, бесконечно малые изменения граничных условий -- приведут к бесконечно малому же изменению решений.

[Фантом]

ewert в сообщении #274207 писал(а):

И, следовательно, бесконечно малые изменения граничных условий -- приведут к бесконечно малому же изменению решений.

Бесконечно малые - да. Но мы же имеем дело, строго говоря, не с дифференциальным уравнением, а с разностной схемой. У которой все погрешности конечны, и соотношение между этими погрешностями из устойчивости (самой по себе) не следует.

[ewert]

Фантом в сообщении #274268 писал(а):

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

Стандартный факт: точность схемы определяется порядком её аппроксимации и её (схемы) устойчивостью.

С последним всё достаточно тривиально, поскольку переменная только одна и поскольку исходный дифференциальный оператор симметричен и знакоопределён (отрицателен). Эти же свойства переносятся и на разностный оператор, если только не изуродовать его специально: раз уж аппроксимация есть, то наименьшее по модулю собственное число (а именно оно определяет норму обратного) примерно совпадает с соотв. собственным числом дифоператора, т.е. с некоторой константой. Следовательно -- задача устойчива.

Это если не принимать во внимание погрешности округления; а если принимать -- то кроме устойчивости надо учитывать ещё и число обусловленности матрицы. Но и с этим всё в порядке: наибольшее по модулю собственное число имеет порядок $h^{-2}$ -- следовательно, такой же порядок имеет и число обусловленности. И при любых практически мыслимых шагах это число не безумно велико. И особенно невелико именно потому, что константа перед второй производной мала: фактически до порядка двадцати узлов число обусловленности будет вообще порядка единицы, и лишь потом начнёт потихонечку возрастать.

[Фантом]

ewert в сообщении #274298 писал(а):

Это если не принимать во внимание погрешности округления; а если принимать -- то кроме устойчивости надо учитывать ещё и число обусловленности матрицы. Но и с этим всё в порядке: наибольшее по модулю собственное число имеет порядок $h^{-2}$ -- следовательно, такой же порядок имеет и число обусловленности. И при любых практически мыслимых шагах это число не безумно велико. И особенно невелико именно потому, что константа перед второй производной мала: фактически до порядка двадцати узлов число обусловленности будет вообще порядка единицы, и лишь потом начнёт потихонечку возрастать.

Логично. Пожалуй, я действительно несколько перестарался.