Clear Filters
Clear Filters

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

66 views (last 30 days)
卓也
卓也 on 10 Dec 2023
Answered: COVAO on 13 Jan 2024
=:角速度 τ:トルク I:慣性行列)の式より得られた角加速度を4次のルンゲクッタを用いて1ステップ時間ごとに数値積分し,角速度,角度を求めたいのですが,実装の方法が分かりません.
Excelファイルの角加速度データを用いて,4次のルンゲクッタの計算方法に関し,,教えていただきたいです.

Answers (1)

COVAO
COVAO on 13 Jan 2024
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

Categories

Find more on 常微分方程式 in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!