Page 127 - 363_
P. 127
128
% z – обчислені значення похідних z(i) =dy(i)/dt.
% h – крок інтегрування
% t – попередній момент часу
% у - попереднє значення вектора змінних стану.
% Вихідні змінні:
% tout – новий момент часу
% yout – обчислене значення вектора змінних стану через крок інтегрування
% Розрахунок проміжних значень похідних
kl = feval(Zpfun, t, у);
k2 = feval(Zpfun, t + h/3, y + h/3 *kl);
k3 = feval(Zpfun, t + 2*h/3, y+h*k2 - h/3*kl);
k4 = feval(Zpfun, t+h, y+h*(k3 + kl - k2));
% Розрахунок нових значень вектора змінних стану
tout = t + h;
yout = у + h*(k l + 3*k2 + 3*k3 + k4)/8;
% Кінець процедури RKO43
Зверніть увагу на такі обставини: звернення до процедури
обчислення правих частин не конкретизовано; ім'я цієї процедури входить до
вхідних змінних процедури інтегрування і конкретизується при зверненні до
останньої; проміжні змінні є векторами-рядками (те саме стосується змінних у,
а також змінних z, що обчислюються у процедурі правих частин.
Процедура правих частин ЗДР маятника. Розглянемо процес
утворення процедури обчислення правих частин ЗДР на прикладі рівняння
маятника, точка підвісу якого поступально переміщується з часом за
гармонічним законом:
J R mgl 1 ( n sin( t )) sin mgl n sin( t ) cos ,
my y mx x
де J – момент інерції маятника; R – коефіцієнт демпфірування; mgl –
опорний маятниковий момент маятника; n my – амплітуда віброперевантаження
точки підвісу маятника у вертикальному напрямку; n mх – амплітуда
віброперевантаження точки підвісу маятника у горизонтальному напрямку;
– кут відхилення маятника від вертикалі; – частота коливань точки підвісу;
х, у – початкові фази коливань (по прискоренню) точки підвісу у
горизонтальному і вертикальному напрямках.
Щоб скласти М-файл процедури обчислення правих частин заданої
системи ЗДР, перш за все слід привести початкову систему ЗДР до форми
Коші. Для цього введемо позначення:
y l = ; y 2 .
Тоді початкове рівняння маятника можна подати у вигляді двох
диференційних рівнянь першого порядку:
