How to normalize unit vector in ode45?

조회 수: 15 (최근 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 CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by