Метод конечных разностей в Mathematica

[grena5]

Добрый день! Следует решить следующую линейную краевую задачу методом конечных разностей в Mathematica:
$y''+2/(x+1)y'+\pi^2y=0$
$y(0)-y'(0)=2-\pi$
$y(1)+y'(1)=-\pi/2-1/4$
Точное решение будет: $(\cos(x\pi)+\sin(x\pi))/(1+x)$
У меня почему-то не правильный ответ получается. Пробовала и заново перебивать и всё несколько дней проверяю, но ошибку найти никак не могу, помогите, пожалуйста найти ошибку. Хоть в какой части программы она.

resh := (p = 2/(x + 1); q = \[Pi]^2; a = 0; A = 2 - \[Pi]; b = 1; 
  B = -\[Pi]/2 - 1/4; n = 5; h = (b - a)/n; f = 0;
  a11 = -1;
  (*Введены начальные данные*)
  (*Разбитие по оси х и вычисление начальных коэффициентов*)
  M = Array[с, {n, 3}];
  For[i = 1, i < n + 1, i++, xi = a + h*i; 
   pi = 2/(xi + 1);
   (*yi+1 + miyi + ki yi-1 = fi*)
   a1 = (2*q*h^2 - 4)/(2 + pi*h); (*mi*)
   a2 = (2 - pi*h)/(2 + pi*h); (*ki*)
   a3 = 3*h^2*f/(2 + pi*h);(*fi*) 
   M[[i]] = {a1, a2, a3}; 
   ]  ;
  (*Print[N[M//MatrixForm,15]];*)
  p1 = 2/(a + h*1 + 1 );
  k1 = (2 - p1*h)/(2 + p1*h);
  m1 = (2*q*h^2 - 4)/(2 + p1*h);
  f1 = 0;
  (*записываеим первое граничное условие, 
  с учетом формулы 1 из блока*)
   (*s,r,...*)
  r = 2 h - 3*a11 + a11*k1;
  s = 4*a11 + a11*m1;
  (*Print[{r,s}];*)
  (*Заключительная формула из блока с учетом второго граничного \
условия*)
  M[[n]] = {a1 - 2*h, a2 + 1, 
    0 - 2*h*B}; (*используя каревое увловие B обнуляем коээфициент \
при yn+1*)
  (*поиск решения системы методом прогонки, приводим к виду ci yi + 
  yi+1 = di*)
  CD = Array[c, {n, 2}];
  CD[[1, 1]] = m1 - s*k1/r;
  CD[[1, 2]] = f1 - (2*h*A + a11*f1)*k1/r;
  For [i = 2, i < n, i++,
   c11 = M[[i, 1]] - M[[i, 2]]/CD[[i - 1, 1]];
   c12 = M[[i, 3]] - CD[[i - 1, 2]]*M[[i, 2]]/CD[[i - 1, 1]];
   CD[[i]] = {c11, c12};
   ];
  (*Print[N[M[[n,3]]-2*h*B-CD[[n-1,2]]*(M[[n,2]]+1)/CD[[n-1,1]],10]];*) 
  v = M[[n, 3]] - 2*h*B - 
    CD[[n - 1, 2]]*(M[[n, 2]] + 1)/CD[[n - 1, 1]];
  u = M[[n, 1]] - 2*h - (M[[n, 2]] + 1)/CD[[n - 1, 1]];
  
  (*Print[CD//MatrixForm]*)
  (*Конец приямой прогонки, начинается вычисление неивестных*)
  Y = Array[c, n + 1];
  Y[[n]] = v/u;
  For[i = n - 1, i > 0, i--,
   Y[[i]] = (CD[[i, 2]] - Y[[i + 1]])/CD[[i, 1]];
   ];
  Y[[n + 1]] = (2*h*A + f1*a11 - s*Y[[1]])/r;
  
  For[i = 2, i < n + 2, i++,
   a1 = Y[[i - 1]];
   OTV[[i]] = a1;
   ];
  OTV[[1]] = Y[[n + 1]];
  
  (*Вывод ответа и проверка на правильность*)
  Print["Ответ:"];
  (*Print[N[Y//MatrixForm,5]] ; *)
  Print[N[OTV // MatrixForm, 5]];
  
  Prov = Array[с, n];
  For [i = 2, i < n, i++,
   pr = Y[[i + 1]] + M[[i, 1]]*Y[[i]] + M[[i, 2]]*Y[[i - 1]] - 
     M[[i, 3]];
   Prov[[i]] = pr;
   ];
  
  Prov[[1]] = Y[[n + 1]]*r + s*Y[[1]] - 2*h*A - f1*a11;
  Prov[[n]] = (M[[n, 2]] + 1)*Y[[n - 1]] + (M[[n, 1]] - 2*h)*Y[[n]] - 
    M[[n, 3]] + 2*h*B;
  (*pr=Y[[n+1]]*r+s*Y[[1]]-2*h*A-f1*a11;*)
  (*pr=(M[[n,2]]+1)*Y[[n-1]]+(M[[n,1]]-2*h)*Y[[n]]-M[[n,3]]+2*h*B;*)
  (*pr=CD[[3,2]]-M[[3,3]]+CD[[2,2]]*M[[3,2]]/CD[[2,1]]; *)
  (*pr=(2*h-3*a11)*Y[[n+1]]+4*a11*Y[[1]]-a11*Y[[2]]-2*h*A; *)
  (*pr=Y[[3]]+M[[2,1]]*Y[[2]]+M[[2,2]]*Y[[1]]-M[[2,3]];*)
  (*Print["Проверка:"];
  Print[N[Prov//MatrixForm,5]]; *)
  
  )

Алгоритм, по которому я писала программу находится в даной методичке

[TOTAL]

Так ошибку найти трудно. Если бы Вы описали здесь алгоритм таким образом, что можно было бы оценить его и оценить то, как он реализован в программе, то шансы найти ошибку заметно увеличились бы.

Кстати, попробуйте краевые условия сразу аппроксимировать на двух точках. Например, при $x=0:$

$$y_0-\frac{1}{1-h} \left(\frac{y_1 - y_0}{h}+\frac{h}{2}\pi^2y_0 \right)=2-\pi$$

[grena5]

Для удобства разберу программу по частям чётко описав, что я делаю:

1) ввожу мои начальные данные, и заменяю производные их приближенными значениями, записывая новые коэффициенты в массив
$y_{i+1} + m_i y_i + k_i y_{i-1} =\varphi _i$
$m_i=\frac{2q_ih^2-4}{2+p_ih}$
$k_i=\frac{2-p_ih}{2+p_ih}$
$\varphi_i=\frac{2h^2f_i}{2+p-ih}$

resh := (p = 2/(x + 1); q = \[Pi]^2; a = 0; A = 2 - \[Pi]; b = 1; 
  B = -\[Pi]/2 - 1/4; n = 5; h = (b - a)/n; f = 0;
  a11 = -1;
  (*Введены начальные данные*)
  (*Разбитие по оси х и вычисление начальных коэффициентов*)
  M = Array[с, {n, 3}];
  For[i = 1, i < n + 1, i++, xi = a + h*i; 
   pi = 2/(xi + 1);
   (*yi+1 + miyi + ki yi-1 = fi*)
   a1 = (2*q*h^2 - 4)/(2 + pi*h); (*mi*)
   a2 = (2 - pi*h)/(2 + pi*h); (*ki*)
   a3 = 3*h^2*f/(2 + pi*h);(*fi*) 
   M[[i]] = {a1, a2, a3}; 
   ]  ;

2) привожу матрицу к трехдиагональному виду используя граничные условия
из первого граничного условия получаю:
$ry_0+sy_1=2hA+\varphi_1 \alpha_1$
$r=2h\alpha_0 - 3 \alpha_1 + k_1 \alpha_1$

$s=4 \alpha_1+ m_1 \alpha_1$

  p1 = 2/(a + h*1 + 1 );
  k1 = (2 - p1*h)/(2 + p1*h);
  m1 = (2*q*h^2 - 4)/(2 + p1*h);
  f1 = 0;
  (*записываеим первое граничное условие, 
  с учетом формулы 1 из блока*)
   (*s,r,...*)
  r = 2 h - 3*a11 + a11*k1;
  s = 4*a11 + a11*m1;

3) Обнуляю коэффициент при $y_{n+1}$ в последнем из системы уравнений (коэффициенты этой системы я на 1ом шаге занесла в массив) с учетом второго граничного условия и получаю в результате:
$\beta_1(k_n+1)y_{n-1}+(m_n\beta _1 -2h \beta_0)y_n=\varphi _n \beta_1 - 2hB$

  M[[n]] = {a1 - 2*h, a2 + 1, 
    0 - 2*h*B}; (

[TOTAL]

grena5 в сообщении #636504 писал(а):

Для удобства разберу программу по частям чётко описав, что я делаю:

1) ввожу мои начальные данные, и заменяю производные их приближенными значениями, записывая новые коэффициенты в массив
$y_{i+1} + m_i y_i + k_i y_{i-1} =\varphi _i$

Что было до этого, как аппроксимировали уравнение?
Как аппроксимировали граничные условия?

-- Сб окт 27, 2012 17:48:50 --

grena5 в сообщении #636504 писал(а):

2) привожу матрицу к трехдиагональному виду используя граничные условия

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

[grena5]

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

[TOTAL]

grena5 в сообщении #636504 писал(а):

1) ввожу мои начальные данные, и заменяю производные их приближенными значениями, записывая новые коэффициенты в массив

Покажите уравнения (не преобразованные!) сразу после замены производных приближенными значениями.

[grena5]

А уравнений у меня нет (в методичке они в явном виде не описаны)
По идее я всё делала четко как в методичке, которую нам дал препод, всё описано на страницах 6-8 вот этой методички https://docs.google.com/open?id=0B8zBQqIOj_AIaFVpWlNxYTlaZjg

[TOTAL]

Дальше только после того, как приведете здесь уравнение, в котором производные аппроксимированы конечными разностями.

[grena5]

Во внутренних узлах сетки производные, входящие в уравнение $y''+p(x)y'+q(x)y=f(x)$ апроксимируются с помощью формул:
$y'(x)=\frac{y_{i+1}-y_{i-1}}{2h}+O(h^2)$
$y'(x)=\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}+O(h^2)$
В точке b мы для апроксимации используем первую приведенную формулу
В точке a применяем формулу:$y'(x_0)= \frac{-y_2+4y_1-3y_0}{2h}+O(h^2)$

[TOTAL]

grena5 в сообщении #636527 писал(а):

Во внутренних узлах сетки производные, входящие в уравнение $y''+p(x)y'+q(x)y=f(x)$ апроксимируются с помощью формул:

Запишите здесь свое уравнение, в котором производные заменены конечными разностями.

[grena5]

Я не понимаю как заменять вторую поризводную

[TOTAL]

grena5 в сообщении #636540 писал(а):

Я не понимаю как заменять вторую поризводную

Вот как: $y''_i=\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}+O(h^2)$

[grena5]

Я сделала это в общем виде и получила четко что в формулах для коэффициентов. Для доказательства этого напишу здесь кусочек моих формул
$y_{i+1}(1/h^2+p/2h)+y_i(-2/h^2+q)+y_{i-1}(1/h^2-p/2h)$
Мне необходимо подставлять сюда значения p=2/(x+1) и $q=\pi^2$

[TOTAL]

grena5 в сообщении #636571 писал(а):

Я сделала это в общем виде и получила четко что в формулах для коэффициентов. Для доказательства этого напишу здесь кусочек моих формул
$y_{i+1}(1/h^2+p/2h)+y_i(-2/h^2+q)+y_{i-1}(1/h^2-p/2h)$
Мне необходимо подставлять сюда значения p=2/(x+1) и $q=\pi^2$

Продолжайте дальше сами, коли не хотите делать то, что я предлагаю.

[grena5]

Спасибо за помощь! Нашла ошибку! Она оказалась совсем не в том месте, где я подозревала