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

[bullvinkle]

Добрый день. Помогите, пожалуйста, исправить программу. Не верно считает коэффициенты матрицы, хотя все вроде должно быть пучком.
Есть задание Решить краевую задачу для обыкновенных дифференциальных уравнений методом конечных разностей.
Краевая задача:
$y''+p(x)y'+q(x)y=f(x)$
$ \left\{ \begin{array}{l} \alpha_0y(a)+\alpha_1y'(a) = A,\\ \beta_0y(b)+\beta_1y'(b)=B, \end{array} \right. $
где
$\alpha_0=5;$ $\alpha_1=0,$
$\beta_0=0;$ $\beta_0=-5;$
$k_1=-1;$ $k_2=-4,8$
$f(x)=k_1-x;$
$q(x)=k_2+x^3 ;$
$p(x)= \sqrt{x^2+k_2^2}$
$A=1 ;$ $B=0.5 ;$
$a=2 ;$ $b=3$
$h = 0.1$ - шаг разбиения, нужен для программы.

Для решения конкретно этих условий, я написал программу, которая находит коэффициенты СЛАУ (и пишет их в D:\1.TXT).
arr[i,j] - матрица; brr[i] - столбец свободных коєфициентов.

const k1=-1;
      k2=-4.8;
      a0=5;
      a1=0;
      b0=0;
      b1=-5;
      AA=1;
      BB=0.5;
var h:real;
    a,b:integer;
    arr:array [0..20,0..20] of real;
    brr:array [0..20] of real;
    i,j:integer;
    x:array [0..20] of real;
    y:array [0..20] of real;
    n:integer;
    fil:text;
function p (var cof:real):real;
 begin
     p:=k1-cof;
 end;
function q (var cof:real):real;
 begin
     q:=k2+(cof*cof*cof);
 end;
function f (var cof:real):real;
 begin
     f:=sqrt(cof*cof+k2*k2);
 end;

BEGIN
     assign(fil,'d:\1.txt');
     rewrite(fil);
     a:=2;
     b:=3;
     h:=0.1;
     n:=(b-a)*10;
     for i:=0 to n do begin
         for j:=0 to n do begin
             arr[i,j]:=0;  end;end;
     arr[0,0]:=a0-3.0*a1/(2.0*h);
     arr[0,1]:=2.0*a1/h;
     arr[0,2]:=a1/(2.0*h);
     brr[0]:=AA;
     for i:=1 to (n-1) do
         begin
           x[i]:=a+i*h;
           arr[i,i-1]:=1.0/(h*h)-p(x[i])/(2.0*h);
           arr[i,i]:=-2.0/(h*h)+q(x[i]);
           arr[i,i+1]:=1.0/(h*h)+p(x[i])/(2.0*h);
           brr[i]:=f(x[i]);
         end;
     arr[n,n-2]:=b1/(2.0*h);
     arr[n,n-1]:=-b1/h;
     arr[n,n]:=b0+3.0*b1/(2.0*h);
     brr[n]:=BB;

    for i:=0 to n do begin
        for j:=0 to n do begin
           write (fil,{'arr[',i,',',j,']:= ',}arr[i,j]:4:2,#09,#09);
        end;
        writeln(fil,#13);
    end;
    writeln(fil,#13);
    writeln(fil,#13);
    writeln(fil,#13);
    for i:=0 to n do begin
          writeln (fil,'brr[',i,']= ',brr[i]:4:2,'  ');
        end;
END.

Писал по вот этой блок-схеме:

В данной блок-схеме: $n$ - количество отрезков разбиения, $a$, $b$ - левая и правая граница отрезка, $\alpha_0, \beta_0, \alpha_1, \beta_1, A, B$ - константы из краевых условий, $h$ - шаг разбиения, $x (0), ... , x (n)$ - система равноудалены точек, $y (0), ..., y (n)$ - приближенные значения решения краевой задачи, $p (x (i)), q (x (i)), f (x ( i))$ - значение функций $p (x), q (x), f (x)$ из дифференциального уравнения в точке $x (i), gauss (a, m, b, y)$ - подпрограмма решения системы линейных алгебраических уравнений методом Гаусса ($a$ - матрица постоянных коэффициентов, $m$ - количество уравнений, $b$ - вектор-столбец свободных членов, $y$ - вектор-столбец неизвестных.
Единственное, что не совпадает, так это нет подпрограммы решения СЛАУ методом Гауса (gaus()). У меня есть отдельная программа, работающая 100% и мне лень переписывать текст каждый раз, как она нужна где-то в др. программе:).

Ответы не совпадают с ответами у преподавателя.
Так же, если у вас, вдруг, есть программа, которая мне поможет, я ее с удовольствием использую. Т.к. преподавателя интересует конкретно ответ, а не то, каким образом я к нему пришел.

[ewert]

Вот эти две формулы явно неверны:

arr[0,2]:=a1/(2.0*h);

arr[n,n-1]:=-b1/h;

Остальное похоже на правду (хотя числовых значений я,естественно, не проверял).

Но, между прочим, ещё пара технических замечаний.

1). Кто ж так пишет:

     a:=2;
     b:=3;
     h:=0.1;
     n:=(b-a)*10;

?
Надо

n:=10;    h:=(b-a)/n;

. И границы интервала задавать как вещественные, а не целые.

2). Решение крайне не рационально. У Вас матрица трёхдиагональна, и решать систему надо методом прогонки, храня в памяти только элементы трёх диагоналей (ну плюс ещё пару элементов из первой и последней строки).

[bullvinkle]

Первую формулу подправил (не углядел минус).
А вот чем Вам не понравилась

arr[n,n-1]:=-b1/h;

?
Вроде все как в схеме.

[ewert]

bullvinkle в сообщении #213927 писал(а):

Вроде все как в схеме.

Ну в схеме и ошибка. Двойка потеряна. Сравните с левым граничным условием -- коэффициенты ж должны быть симметричны.

[bullvinkle]

Спасибо.
Т.е. нужно

arr[n,n-1]:=-b1/(2.0*h);

Завтра пойду сдавать. Потом отвечу верно ли.

[ewert]

bullvinkle в сообщении #213934 писал(а):

Спасибо.
Т.е. нужно

arr[n,n-1]:=-b1/(2.0*h);

Да наоборот, конечно -- вверх двойку. Какие коэффициенты у трёхточечной производной?

[bullvinkle]

Двойку отправил в числитель.
Термин трехточечная производная вижу впервые, не знаю как ответить.
С учетом всех исправлений у меня вот такой результат:
http://www.sharemania.ru/0144475 , если Вы это имели ввиду.

[ewert]

bullvinkle в сообщении #213967 писал(а):

Термин трехточечная производная вижу впервые, не знаю как ответить.

"Эн-точечные разностные производные" -- это приближённые производные, вычисляемые по эн узлам. Возможно, это и жаргон, но вполне общеупотребительный.

$$y'(x_0)={-3y_0+4y_1-y_2\over 2h}+O(h^2);$$

$$y'(x_n)={y_{n-2}-4y_{n-1}+3y_n\over 2h}+O(h^2).$$

Файл у меня: 1) не загружается; 2) килобайт спокойно можно было бы и непосредственно в сообщение втиснуть; 3) результаты мне всё равно ничего не скажут.

[bullvinkle]

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

Цитата:

2) килобайт спокойно можно было бы и непосредственно в сообщение втиснуть;

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

[ewert]

Ну у Вас же текстовой файл. Копипастните его и тисните в окошко ответа. Желательно в теге "код", чтобы не нарушить форматирование.

А коэффициенты -- они и есть коэффициенты, те, что в числителях, скажем.

[bullvinkle]

5.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      

115.50   -195.54      84.50      0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      

0.00      116.00      -194.15      84.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      

0.00      0.00      116.50      -192.63      83.50      0.00      0.00      0.00      0.00      0.00      0.00      

0.00      0.00      0.00      117.00      -190.98      83.00      0.00      0.00      0.00      0.00      0.00      

0.00      0.00      0.00      0.00      117.50      -189.17      82.50      0.00      0.00      0.00      0.00      

0.00      0.00      0.00      0.00      0.00      118.00      -187.22      82.00      0.00      0.00      0.00      

0.00      0.00      0.00      0.00      0.00      0.00      118.50      -185.12      81.50      0.00      0.00      

0.00      0.00      0.00      0.00      0.00      0.00      0.00      119.00      -182.85      81.00      0.00      

0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      119.50      -180.41      80.50      

0.00      0.00      0.00      0.00      0.00      0.00      0.00      0.00      -25.00      100.00      -75.00      

brr[0]= 1.00  
brr[1]= 5.24  
brr[2]= 5.28  
brr[3]= 5.32  
brr[4]= 5.37  
brr[5]= 5.41  
brr[6]= 5.46  
brr[7]= 5.51  
brr[8]= 5.56  
brr[9]= 5.61  
brr[10]= 0.50  

Вот текстовый файл. Матрица и столбец свободных коэффициентов.

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

[ewert]

bullvinkle в сообщении #214035 писал(а):

Вот текстовый файл. Матрица и столбец свободных коэффициентов.

Ну и что? Матрица как матрица, такой её структура и должна быть. А конкретные числа мне, естественно, ни о чём не говорят.

bullvinkle в сообщении #214035 писал(а):

, но не знаю где и как мне, в конкретно моей задаче, найти эти приближённые производные, вычисляемые по эн узлам

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

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

[bullvinkle]

Только вчера сдал эту лабу (преподаватель где-то затерялся).
Все правильно, ответы совпадают.
Большое спасибо за инструктаж=).

[21rus]

объясните откуда берется cof "p:=k1-cof;" здесь и так далее? где оно инициализируется?ни как не могу понять .. вроде по формулам там идет X

[EtCetera]

В данном случае $p$ - функция, а не переменная. А потому инициализироваться не должна...