Практикум по методам решения задачи Коши для систем ОДУ. Даутов P.З. - 18 стр.

18                                                                   Задача Коши для системы ОДУ


   Cхематически заголовок программы RKq можно определить в
виде
             f unction [T, Y ] = RKq(f, ab, ua, N ).
Здесь в квадратных скобках указываются выходные параметры, а в
круглых — входные.1)
     Для заданной матрицы Y под Y (i, :) будем понимать i-тую строку
Y , а для заданных одномерных массивов x, y, длины n, и чисел c, d,
под x = y, z = c ∗ x + d ∗ y, y = 0 будем понимать реализацию
покомпонентных операций
               xi = yi ,       z i = c ∗ xi + d ∗ y i ,         yi = 0,        i = 1, 2, . . . , n.
Будем считать также, что команда Y = zeros(m, n) определяет мат-
рицу размера m × n, состоящую из одних нулей; под x(i) будем по-
нимать i-тый элемент массива (вектора строки или столбца). Ниже
приведена в схематической форме реализация 3-х стадийного метода
(для примера); числа типа alpha2, beta21 и т.д. должны быть заме-
нены на конкретные значения.
f u n c t i o n [ T,Y]=RK3( f , ab , ua ,N)
T=z e r o s ( 1 ,N+1);
Y=z e r o s (N+1,n ) ;

t=ab ( 1 ) ;
y=ua ;
h=(ab (2) − ab ( 1 ) ) /N;

T(1)= t ;
Y( 1 , : ) = y ;
f o r i =1:N
    k1=f ( t , y ) ;                 % i n MatLab : k1=f e v a l ( f , t , y ) ;

      yp=y+( b e t a 2 1 *h) * k1 ;
      k2=f ( t+a l p h a 2 *h , yp ) ; % i n MatLab : k2=f e v a l ( f , t+a l p h a 2 *h , yp ) ;

      yp=y+( b e t a 3 1 *h) * k1+( b e t a 3 2 *h) * k2 ;
      k3=f ( t+a l p h a 3 *h , yp ) ; % i n MatLab : k3=f e v a l ( f , t+a l p h a 3 *h , yp ) ;

     y=y+h * ( p1 * k1+p2 * k2+p3 * k3 ) ;

  t=t+h ;
  T( i +1)=t ;
  Y( i +1 ,:)= y ;
end
   1)
        такой синтаксис языка программирования принят в MatLab. Мы будем его придерживать-
ся.