Integrate 61 x 3 tabular data using Runge Kutta

I have a 61 x 3 velocity data and wanted to integrate using Rung Kutta integrator? I know how to integrate when I have a function but I am a bit confused when it is tabular data. Integration step is . Data is attached. Any help is apperciated. Thanks

Answer 1

[Ameer Hamza]

If you are okay with using the built-in function, then use ode45, which is an implementation of the RK (4,5) algorithm.

data = load('angular.velocity.data.txt');

v1 = data(:,1); % First column
v2 = data(:,2); % Second column
v3 = data(:,3); % Third column
time = 0:0.016:0.016*(numel(v1)-1);
IC = 0; % initial displacement is 0

[~,x1] = ode45(@(t,x) odeFun(t, x, time, v1), time, 0);
[~,x2] = ode45(@(t,x) odeFun(t, x, time, v2), time, 0);
[~,x3] = ode45(@(t,x) odeFun(t, x, time, v3), time, 0);

%% plotting
subplot(1,3,1)
plot(t,v1,t,x1);
legend({'Velocity', 'displacement'})
subplot(1,3,2)
plot(t,v2,t,x2);
legend({'Velocity', 'displacement'})
subplot(1,3,3)
plot(t,v3,t,x3);
legend({'Velocity', 'displacement'})

function v = odeFun(t, x, time, vx)
%     v = interp1(time, vx, t, 'linear'); % linear interpolation
%     v = interp1(time, vx, t, 'pchip'); % pchip interpolation
    v = interp1(time, vx, t, 'makima'); % makima interpolation
end

[HN]

Thank you Ameer Hamza,

This method also works to integrate velocity inside the loop before all data's is stored ? I used for loop to compute v1, v2 and v3.

Thanks a lot

[Ameer Hamza]

I am glad to be of help.

Yes, you can also use a for-loop to integrate the velocity values.

[HN]

%I've tried to use the program inside the for loop and it gives me the following error. Thanks agian 
Error using odearguments (line 21)
When the first argument to ode45 is a function handle, the tspan argument must have at
least two elements.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in MainPRS20200405 (line 251)
    [~,x1] = ode45(@(t,x) odeFun(t, x, time, vvx), time, 0);

[Ameer Hamza]

Can you show your code?

[HN]

   %vector v runs inside the for loop which is indexed by k
   t=1:delt:1;
   for k=1:length(t)
       v(:,k)=[4*cos(2*pi*t(k));4*sin(2*pi*t(k));10*cos(2*pi*t(k))]*0;%Jv*dgc;
        vx(k)=v(1,k);
        vy(k)=v(2,k);
        vz(k)=v(3,k);
    
        time = 0:0.016:0.016*(numel(vx)-1);
        IC = 0; % initial displacement is 0
        [~,x1] = ode45(@(t,x) odeFun(t, x, time, vx), time, 0);
        [~,x2] = ode45(@(t,x) odeFun(t, x, time, vy), time, 0);
        [~,x3] = ode45(@(t,x) odeFun(t, x, time, vz), time, 0);
   end

[Ameer Hamza]

No, this is not the correct approach. Why do you want to write a for-loop like this? What is the goal here?

[HN]

Actually the code here with me is much complicated than the above. So, I sent you simpler one to show what I am doing . Here is the acutal problem I wanted to solve. Matrix G is made of . So must be integrated to run the next loop and update G. Thats what I want.

Thanks you

   t=1:0.016:1;
for k=1:length(t)
    .
    .
    .
    G=inv(Jq)*Gx;   
    detR(k)=det([[G(1:3,:)];[G(4:6,:)]]'); % Determinant   
    
       
    v(:,k)=[4*cos(2*pi*t(k));4*sin(2*pi*t(k));10*cos(2*pi*t(k))]*0;%Jv*dgc;           

    w(:,k)=[-sin(2*pi*t(k)); cos(2*pi*t(k));sin(2*pi*t(k))]*10;
    
    st=[vx(k);vy(k);vz(k);wx(k);wy(k);wz(k)];
    MotionModifierM=[eye(6)-transpose(G(7:9,:))*pinv((transpose(G(7:9,:))))];
    sta=MotionModifierM*st;%([eye(6)-transpose(G(7:9,:))*pinv((transpose(G(7:9,:))))]*st(:,k)); % modfied velocity

    dq(:,k)=G*sta; 
    
    dd1(k)=dq(1,k);
    dd2(k)=dq(2,k);
    dd3(k)=dq(3,k);
    
    time = 0:0.016:0.016*(numel(vx)-1);
    IC = 0; % initial displacement is 0
    [~,x1] = ode45(@(t,x) odeFun(t, x, time, dd1), time, 0);
    [~,x2] = ode45(@(t,x) odeFun(t, x, time, dd2), time, 0);
    [~,x3] = ode45(@(t,x) odeFun(t, x, time, dd3), time, 0);
end

[Ameer Hamza]

You cannot use ode45 like that. The code I have written in my answer requires you to have all the values in the elements in a vector. Why don't you first calculate all the values of v and w and then apply ode45.

[HN]

I used your code to integrate v and w. But I still have to integrate dq(:,k) that should be used to update G in each iteration. Thats why I cameup with the new question. I integrate dq using Euler's integrartion but integration of v& w with Rung Kutta and integration of dq with Eulers should be compared to check the error . However, I felt it is not possible to get proper error with diffeent integrtion method.

[Ameer Hamza]

Can you show you actual equations in mathematical form?

[Ameer Hamza]

Yes, you can share that too.

[Ameer Hamza]

Value of which variables change with the value of integration? See this example: https://www.mathworks.com/help/releases/R2020a/matlab/ref/ode45.html#bu3l43b to see how to handle parameters which vary with the input variable.

[HN]

Ameer Hamza ,

I need your help in similr problem,

S  % 3 by 100 data  
W=skew(S) % 4 by 4 marix 
dq=0.5*(W*q')';  % have only initial q 

I wanted to integrate dq to get q. Can you help me ?

Answer 2

[darova]

• YOu can create function using interp1 or spline
• You can write your own solver

Simple solver Euler method

for i = 1:length(velocity)
    position(i+1) = position(i) + dt*velocity(i);
end

[HN]

Thank you!