Метод конечных разностей 4 порядка, для второй производной..

[Ruski]

Здравствуйте!
Очень нужен алгоритм метода конечных разностей 4- го порядка, для решения краевой задачи
вида:
$\ddot{y}-p(x)y=f(x)$

$x \in [a, b] $
$p(x) \ge 0$
$y(a)=A$
$y(b)=B$

нужен метод 4-го порядка...
Я почитал информацию на сайте: http://ad.cctpu.edu.ru/Math_method/math/12.htm
Но там для 2-го порядка(
Так же в учебнике по численным методам Н.С. Бахвалова прочитал, что для 4-го порядка
справедлива схема :
$h = (b-a)/N$
N - число узлов

$\frac{y_{i+1} -2y_{i} +y_{i-1} }{h^{2} } -p_{i} y_{i} -\frac{1}{12} ((p_{i+1} y_{i+1} +f_{i+1} )-2(p_{i} y_{i} +f_{i} )+(p_{i-1} y_{i-1} +f_{i-1} ))=f_{i} $


Аналогично алгоритму, описанному на сайте (http://ad.cctpu.edu.ru/Math_method/math/12.htm)
я попытался реализовать метод 4-го порядка...

Получил алгоритм:


Система выглядит так:
$\noindent \[m_{i} =-\left(\frac{10p_{i} +\frac{24}{h^{2} } }{p_{i+1} -\frac{12}{h^{2} } } \right)\begin{array}{cc} {} & {} \end{array};\begin{array}{cc} {} & {} \end{array}ni=\frac{\left(\frac{12}{h^{2} } -p_{i-1} \right)}{\left(\frac{12}{h^{2} } -p_{i+1} \right)} \] $

$\[y_{i+1} +y_{i} m_{i} +y_{i-1} n_{i} =(f_{i+1} +10f_{i} +f_{i-1} )\frac{h^{2} }{12-h^{2} p_{i+1} } \] $

$ a_{0} y(x_{0} )=A $
$ b_{0} y(x_{n} )=B $

1) ввожу N, a, b, a0, b0, A, B

2) определяю $h = (b-a)/N$

3) определяю сетку $x_i$ с шагом h
где
$x_0 = a$
$x_n = b$

4) определяю коэффициенты $c_i$ и $d_i$
с помощью полученных соотношений:
$c_{1} =\frac{1}{m_{1} } $

$d_{1} =(f_{2} +10f_{1} +f_{0} )-n_{1} \frac{A}{a_{0} } $

$c_{i} =\frac{1}{m_{i} -n_{i} c_{i-1} } $

$d_{i} =(f_{i+1} +10f_{i} +f_{i-1} )-n_{i} c_{i-1} d_{i-1}$

$i=2,3,...,n$

5) определив все коэффициенты $c_i$ и $d_i$
нахожу все $y(x_i)$

$y_{n} =\frac{B}{b_{0} } $

$y_{0} =\frac{A}{a_{0} } $

$y_{i} =c_{i} (d_{i} -y_{i+1} )$

$i=n-1,\begin{array}{cc} {} & {n-2,\begin{array}{cc} {...,} & {1} \end{array}} \end{array} $

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

[ewert]

Ruski в сообщении #579225 писал(а):

Аналогично алгоритму, описанному на сайте (http://ad.cctpu.edu.ru/Math_method/math/12.htm)
я попытался реализовать метод 4-го порядка...

Там какое-то безумие -- очень длинно и совершенно не приходя в сознание.

Схема из Бахвалова, которую Вы привели -- действительно 4-го порядка (это т.наз. метод Нумерова). Подобные приведены, кажется, тоже верно. А вот дальше начинается какое-то совершенно ненужное колдовство.

Надо чётко понимать, что составление системы уравнений и её решение -- это независимые подзадачи, и программировать их следует независимо. Метод прогонки, о котором с таким придыханием говорят http://ad.cctpu.edu.ru/Math_method/math/12.htm -- это в данном случае всего-навсего стандартный метод Гаусса, применённый к трёхдиагональной системе:

$$\begin{pmatrix}d_1&e_1&0&0&\cdots&\vrule&b_1 \\ c_2&d_2&e_2&0&\cdots&\vrule&b_2 \\ 0&c_3&d_3&e_3&\cdots&\vrule&b_3 \\ 0&0&c_4&d_4&\cdots&\vrule&b_4 \\ \vdots&\vdots&\vdots&\vdots&\ddots&\vrule&\vdots&\end{pmatrix} \rightarrow \begin{pmatrix}1&a_1&0&0&\cdots&\vrule&\tilde b_1 \\ 0&1&a_2&0&\cdots&\vrule&\tilde b_2 \\ 0&0&1&a_3&\cdots&\vrule&\tilde b_3 \\ 0&0&0&1&\cdots&\vrule&\tilde b_4 \\ \vdots&\vdots&\vdots&\vdots&\ddots&\vrule&\vdots&\end{pmatrix} \rightarrow \begin{pmatrix}1&0&0&0&\cdots&\vrule&y_1 \\ 0&1&0&0&\cdots&\vrule&y_2 \\ 0&0&1&0&\cdots&\vrule&y_3 \\ 0&0&0&1&\cdots&\vrule&y_4 \\ \vdots&\vdots&\vdots&\vdots&\ddots&\vrule&\vdots&\end{pmatrix}.$$

Вот и программируйте это как отдельную процедуру с четырьмя входными массивами: $\vec e,\vec d,\vec c,\vec b$, одним выходным $\vec y$ и внутренним рабочим $\vec a$. Из-за трёхдиагональности рабочие формулы просты и достаточно очевидно выписываются; готовые формулы см., скажем, в http://e-lib.gasu.ru/eposobia/metody/R_1_3.html. Там, правда, немножко другие обозначения: иксы вместо игреков, в качестве внутреннего рабочего массива $\tilde b$ используется входной $b$ и ещё кое-что; но значения это не имеет.

И ещё: надо помнить, что разностные уравнения -- только для внутренних узлов: $y_1,y_2,\ldots,y_{n-1}$. Значения $y_1$ и $y_n$ заданы по условию, и содержащие их слагаемые в первом и последнем уравнении надо перед решением системы перенести в правую часть, т.е. вычесть из соотв. $f$.

[Ruski]

Спасибо большое Вам за ответ.)
Но я не совсем пойму, что содержат в себе входные массивы $\vec e, \vec d, \vec c, \vec b $
применительно
к системе:
$\noindent \[m_{i} =-\left(\frac{10p_{i} +\frac{24}{h^{2} } }{p_{i+1} -\frac{12}{h^{2} } } \right)\begin{array}{cc} {} & {} \end{array};\begin{array}{cc} {} & {} \end{array}ni=\frac{\left(\frac{12}{h^{2} } -p_{i-1} \right)}{\left(\frac{12}{h^{2} } -p_{i+1} \right)} \] $

$\[y_{i+1} +y_{i} m_{i} +y_{i-1} n_{i} =(f_{i+1} +10f_{i} +f_{i-1} )\frac{h^{2} }{12-h^{2} p_{i+1} } \] $

$ a_{0} y(x_{0} )=A $
$ b_{0} y(x_{n} )=B $

[ewert]

Ruski в сообщении #579284 писал(а):

$\[y_{i+1} +y_{i} m_{i} +y_{i-1} n_{i} =(f_{i+1} +10f_{i} +f_{i-1} )\frac{h^{2} }{12-h^{2} p_{i+1} } \] $

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

[Ruski]

Спасибо Вам ewert!!! Все заработало! ошибка была в формуле для $m_i$, надо брать все выражение с плюсом! и все ОК!!!

[ewert]

(Оффтоп)

Ruski в сообщении #579408 писал(а):

надо брать все выражение с плюсом!

Хм. А я этого даже и не заметил. Действительно, знак там был явно перепутан. Так что вся моя "помощь" оказалась впустую -- всего лишь спровоцировал встрепыхнуться...