4次のルンゲクッタ法に関して

= （：角速度　τ：トルク　I:慣性行列）の式より得られた角加速度を4次のルンゲクッタを用いて1ステップ時間ごとに数値積分し，角速度，角度を求めたいのですが，実装の方法が分かりません．

Excelファイルの角加速度データを用いて，4次のルンゲクッタの計算方法に関し,，教えていただきたいです．

Answer 1

[covao]

MATLABのodeソルバーの関数で一般的に使われるode45関数を用いて各加速度から角速度、角度を計算する例です。

(スクリプト作成に生成AIを利用しています)

角加速度データは、サンプリング時間0.01[sec]、単位:[rad/sec^2]を仮定しています。

% Load angular acceleration data from Excel file
excelFile = 'https://jp.mathworks.com/matlabcentral/answers/uploaded_files/1564584/Book1.xlsx'; % Set the path of the Excel file
angularAccelData = readmatrix(excelFile);

% Generate time data
dt = 0.01; % Assume time step
timeData = (0:dt:dt*(length(angularAccelData)-1))';

% Define angular acceleration as a function of time
angularAccelFunc = @(t) interp1(timeData, angularAccelData, t, 'linear', 'extrap');

% Initial conditions
initialConditions = [0; 0]; % Initial angle 0 radians, initial angular velocity 0 rad/sec

% Solve the ordinary differential equation using ode45 solver
[t, y] = ode45(@(t, y) angularMotionODE(t, y, angularAccelFunc), timeData, initialConditions);
theta = y(:, 1);
omega = y(:, 2);
alpha = angularAccelFunc(t);

% Plotting the results
% Alpha (Angular Acceleration)
subplot(3,1,1);
plot(t, alpha);
xlabel('Time [sec]');
ylabel('Alpha [rad/sec^2]');

% Omega (Angular Velocity)
subplot(3,1,2);
plot(t, omega); 
xlabel('Time [sec]');
ylabel('Omega [rad/sec]');

% Theta (Angle)
subplot(3,1,3);
plot(t, theta);
xlabel('Time [sec]');
ylabel('Theta [rad]');

% Function to solve the ordinary differential equation
function dy = angularMotionODE(t, y, angularAccelFunc)
    dy = zeros(2,1);
    dy(1) = y(2); % y(1) is angle, y(2) is angular velocity
    dy(2) = angularAccelFunc(t); % Angular acceleration
end