How to normalize unit vector in ode45?

Question

Sam H 2018-5-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/402069-how-to-normalize-unit-vector-in-ode45

评论： 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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-5-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/402069-how-to-normalize-unit-vector-in-ode45#answer_321519

在 MATLAB Online 中打开

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 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 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