How to normalize unit vector in ode45?

1 回表示 (過去 30 日間)
Sam H
Sam H 2018 年 5 月 23 日
コメント済み: Torsten 2018 年 12 月 7 日
for example I want to solve a unit vector rotating around an axis. the kinematics would be: i_dot = w cross i, where 'cross' means cross product, ' i' is the unit vector and ' w' is angular velocity vector. Then ode45 is used to solve this equation. However, because of the numerical integration, length of the unit vector will drift away from unity, how can I normalize the vector in every loop? Thanks! Example code:
clear;clc
global w
w = [0;0;1]; % rotates around z axis, with angular velocity 1 rad/s
i0 = [1/2;1/2;1/sqrt(2)]; % initial unit vector
tspan = [0,100]; % time span
[t,i] = ode45(@exampleFun,tspan,i0); % use ode45 to solve for new unit vector
function idot = exampleFun(t,i)
global w
idot = cross(w,i);
end

採用された回答

Torsten
Torsten 2018 年 5 月 23 日
function main
w = [0;0;1]; % rotates around z axis, with angular velocity 1 rad/s
i0 = [1/2;1/2;1/sqrt(2)]; % initial unit vector
tspan = [0,100]; % time span
opts = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,i] = ode45(@(t,i)exampleFun(t,i,w),tspan,i0,opts); % use ode45 to solve for new unit vector
absi=i(:,1).^2+i(:,2).^2+i(:,3).^2;
plot(t,absi)
end
function idot = exampleFun(t,i,w)
idot = cross(w,i);
end
Better ?
Best wishes
Torsten.
  2 件のコメント
roya afshar
roya afshar 2018 年 12 月 7 日
if idot = cross(w,i);is a second order differential equation. which changes do we need in these codes?
D2i=cross(w,i)
Torsten
Torsten 2018 年 12 月 7 日
function idot = exampleFun(t,i,w)
D2i = cross(w,i);
idot = zeros(6,1);
idot(1) = i(2);
idot(2) = D2i(1);
idot(3) = i(4);
idot(4) = D2i(2);
idot(5) = i(6);
idot(6) = D2i(3);
end

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDownloads についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by