Решение одномерной краевой задачи методом разностных схем

[fretyno]

Добрый день. Решение есть, но оно где-то неверное, не могу понять где. Необходимо решить задачу на отрезке [0, 1]. Шаг сетки h = 0.01. Количество разбиений N = 1/h = 100.
Итого мы имеем 101 узел сетки и 101 уравнение.
$ $$\left\{ \begin{array}{rcl} -4\frac{d^2u}{dx^2}+\sin x * u(x)=\sin (\frac{\pi}{2}x) \\ \frac{du}{dx}(0)=0 \\ u(1)=2$ \\ \end{array} \right.$$$

Первым шагом заменяю производную 2-го порядка разностным отношением:
$\frac{d^2u}{dx^2}$\approx\frac{u(x+h)-2u(x)+(x-h)}{h^2}$$
Также заменяю производную 1-го порядка на начальном условии разностным отношением вперед и выполняю преобразования над ним:
$\frac{du}{dx}$\approx\frac{u(x+h)-u(x)}{h}$$
$\frac{u(x+h)-u(x)}{h}=0$
$u(x)=u(x+h)$ - т.к. это граничное условие для x[0] узла, значит мы имеем решение для u[0] и u[1] (1-го и 2-го уравнения) в узлах x[0] и x[1] соответственно.
Формула нахождения решения для u(0) выглядит следующим образом, отсюда выразим $\chi[0] и \nu[0]$
${u}_{0}={\chi}_{0}*{u}_{1}+{\nu}_{0}$ - отсюда имеем:
{\chi}_{0}=1, {\nu}_{0}=0


После всех преобразований имеем такую систему:
$ $$\left\{ \begin{array}{rcl} -4{u}_{j+1}+{u}_{j}*(8+{h}^{2}\sin x) - 4{u}_{j-1}={h}^{2}\sin(\frac{\pi}{2}x) \\ u_0=u_1=0 \\ u_N=2$ \\ \end{array} \right.$$$

Далее методом прогонки я решаю СЛАУ из 98 уравнений, т.к. решения u(0), u(1) и u(n) известны. Решаю с помощью C++(если нужно, код прилагаю), но результаты мне не очень нравятся. Т.к. значение u(n) не равно 2.
Пожалуйста, подскажите, в чем ошибка.

Используется синтаксис C++

#include <iostream>
#include <cmath>
using namespace std;
int main() {
    double h = 0.01;
    int N = 100; int k=0;
    double x[100] = {0}, f[100] = {0}, u[100] = {0};
    double X[100] = {0}, V[100] = {0};
    double C[100] = {0};
    double A = -4.0, B = -4.0;
    u[N] = 2;
    u[0] = u[1] = 0;
    X[0] = 1;
    V[0] = 0;
    cout << "du/dx(0) =  " << u[0] << ",  u[N] = " << u[N] <<", h = " << h << endl << endl;
 
    for ( int i = 0; i <= N; i++) {
        x[i] = i*h;
        f[i] = h*h*(sin(((M_PI)/2)*x[i]));
        C[i] = 8+h*h*sin(x[i]);
    }
 
    for (int i = 0; i < N; i++) {
        if (abs(C[i]) < abs(A)+abs(B))
            cout << "Uslovie sxodimosti ne vipolneno!" << endl;
 
    }
 
    for (int i = 2; i <= N; i++) {
        X[i] = B/(C[i]-A*X[i-1]);
        V[i] = (A*V[i-1]-f[i])/(C[i]-A*X[i-1]);
        X[N] = A/C[N];
        V[N] = -f[N]/C[N];
    }
 
    for ( int i = N; i >= 2; i--) {
        u[i] = X[i]*u[i+1]+V[i];
       }
 
    cout << "i | " << "    X    | " << "    V    | " <<  "     u | " << endl;
    cout << "----------------------------------" << endl;
 
    for (int i = 0; i < N; i++)
        cout << i << " | " << X[i] << " | " << V[i] << " | " << u[i] << " | " << endl;
 
    cout << N << " | " << X[N] << " | " << V[N] << " | " << u[N] << " | " << endl;
    cout<<"Proverka: "<<endl;
    cout<< A*u[49]-C[50]*u[50]+B*u[51] << " = " << f[50] <<endl;
}

[Karan]

Наберите формулы местным редактором, поправьте код и объясните, чем именно Вам не нравится результат.

[Karan]

Тема перемещена из форума «Математика (общие вопросы)» в форум «Карантин»
по следующим причинам:

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы);
- неправильно набран код

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

[Karan]

Тема перемещена из форума «Карантин» в форум «Computer Science»

-- 20.11.2016, 13:48 --

Во-первых, у Вас массивы длины 100, $N = 100$, поэтому когда у Вас в цикле $i = N$, у Вас выход за границы массива.
Во-вторых, я не вижу, где у Вас в коде кусок, отвечающий за граничные условия. Я так понимаю, где-то должно устанавливаться f[N] = 2 или u[N] = 2

[fretyno]

Karan в сообщении #1170316 писал(а):

Тема перемещена из форума «Карантин» в форум «Computer Science»

-- 20.11.2016, 13:48 --

Во-первых, у Вас массивы длины 100, $N = 100$, поэтому когда у Вас в цикле $i = N$, у Вас выход за границы массива.
Во-вторых, я не вижу, где у Вас в коде кусок, отвечающий за граничные условия. Я так понимаю, где-то должно устанавливаться f[N] = 2 или u[N] = 2

Ну так индексы массивов с 0 начинаются, а не с 1, если в нем последний индекс 100, то включая нулевой в нем 101 элемент.
Оно и устанавливается там так.

[Xaositect]

fretyno в сообщении #1170319 писал(а):

Ну так индексы массивов с 0 начинаются, а не с 1, если в нем последний индекс 100, то включая нулевой в нем 101 элемент.

Если массив объявлен как double x[100], то в нем последний индекс 99.
А вот в последнем цикле вообще идет обращение к u[101] в первой итерации:

    for ( int i = N; i >= 2; i--) {
        u[i] = X[i]*u[i+1]+V[i];
       }

[fretyno]

Xaositect в сообщении #1170322 писал(а):

fretyno в сообщении #1170319 писал(а):

Ну так индексы массивов с 0 начинаются, а не с 1, если в нем последний индекс 100, то включая нулевой в нем 101 элемент.

Если массив объявлен как double x[100], то в нем последний индекс 99.
А вот в последнем цикле вообще идет обращение к u[101] в первой итерации:

    for ( int i = N; i >= 2; i--) {
        u[i] = X[i]*u[i+1]+V[i];
       }

Тогда почему у меня не возникает ошибки выхода за границу массива?

[Xaositect]

fretyno в сообщении #1170328 писал(а):

Тогда почему у меня не возникает ошибки выхода за границу массива?

Ее нет в большинстве реализаций C++. Он просто тихо берет значение из соседней области памяти, там находится либо мусор, либо соседний массив.

[fretyno]

Xaositect в сообщении #1170330 писал(а):

fretyno в сообщении #1170328 писал(а):

Тогда почему у меня не возникает ошибки выхода за границу массива?

Ее нет в большинстве реализаций C++. Он просто тихо берет значение из соседней области памяти, там находится либо мусор, либо соседний массив.

Ладно, программа меня не сильно сама интересует. Важен ход решения до кода. Верно ли я все сделал?

[Pphantom]

fretyno в сообщении #1170223 писал(а):

После всех преобразований имеем такую систему:

$$\left\{ \begin{array}{rcl} -4{u}_{j+1}+{u}_{j}*(8+{h}^{2}\sin x) - 4{u}_{j-1}={h}^{2}\sin(\frac{\pi}{2}x) \\ u_0=u_1=0 \\ u_N=2$ \\ \end{array} \right.$$

Далее методом прогонки я решаю СЛАУ из 98 уравнений, т.к. решения u(0), u(1) и u(n) известны.

Вы каким-то хитрым образом сделали вывод, что $u_0=u_1=0$, который ниоткуда не следует. Из нулевого значения производной на левой границе можно получить только то, что $u_0=u_1$. В том же виде, к которому Вы переделали задачу, она оказывается задачей Коши с условиями, заданными на левой границе, которые, естественно, могут оказаться согласованными с условием на правой границе только в том случае, если Вам очень сильно повезет.

[fretyno]

Pphantom в сообщении #1170858 писал(а):

Вы каким-то хитрым образом сделали вывод, что $u_0=u_1=0$, который ниоткуда не следует. Из нулевого значения производной на левой границе можно получить только то, что $u_0=u_1$. В том же виде, к которому Вы переделали задачу, она оказывается задачей Коши с условиями, заданными на левой границе, которые, естественно, могут оказаться согласованными с условием на правой границе только в том случае, если Вам очень сильно повезет.

Пересмотрел размеры массивов и индексы. Убрал что u[0]=u[1]=0. Вроде теперь правильно. Для наглядности построил график по значениям, думаю, что-то такое и должно было выйти.

#include <iostream>
#include <cmath>
using namespace std;
int main() {
    double h = 0.01;
    int N = 100;
    double x[102] = {0}, f[102] = {0}, u[102] = {0};
    double X[102] = {0}, V[102] = {0};
    double C[102] = {0};
    double A = -4.0, B = -4.0;
    u[N] = 2;
    X[0] = 1;
    V[0] = 0;
    cout << "du/dx(0) =  " << u[0] << ",  u[N] = " << u[N] <<", h = " << h << endl << endl;

    for (int i = 0; i <= N; i++) {
        x[i] = i*h;
        f[i] = h*h*(sin(((M_PI)/2)*x[i]));
        C[i] = 8+h*h*sin(x[i]);
    }

    for (int i = 0; i <= N; i++)
        if (abs(C[i]) < abs(A)+abs(B))
            cout << "Uslovie sxodimosti ne vipolneno!" << endl;

    for (int i = 1; i <= N; i++) {
        X[i] = B/(C[i]-A*X[i-1]);
        V[i] = (A*V[i-1]-f[i])/(C[i]-A*X[i-1]);
    }

    for (int i = N-1; i >= 0; i--)
        u[i] = X[i]*u[i+1]+V[i];

    cout << "i | " << "    X    | " << "       V      | " <<  "   u   | " << endl;
    cout << "----------------------------------" << endl;
    cout << 0 << " | " << X[0] << "         | " << V[0] << "            | " << u[0] << " | " << endl;
    for (int i = 1; i < N; i++)
        cout << i << " | " << X[i] << " | " << V[i] << " | " << u[i] << " | " << endl;

    cout << N << " | " << X[N] << " | " << V[N] << " | " << u[N] << " | " << endl;
    cout<<"Proverka: "<<endl;
    cout<< A*u[49]-C[50]*u[50]+B*u[51] << " = " << f[50] <<endl;

}

[Pphantom]

Нет, это явно неправильный результат, ищите ошибку.

[fretyno]

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

[Pphantom]

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

[fretyno]

Pphantom в сообщении #1171600 писал(а):

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