Equation of Motion Solution via Runge Kutta

Hi !

This is my first attemp to use MatLab. For my project, I need to plot a q (input angle) vs. t(time) graph. The fi equations are for linkers' angles and g and g' functions are for velocity and acceleration influence coefficients. The problem about the runge kutta solution is that, my angles, influence coefficient, the force that is shown with Fc, and the input Torque (T) are changing with the q and qdot. So I try to change them for every q and qdot in the for loop. Also since my Q,C and J fuctions are changing, the q double dot is effected , again I put that q double that function in for loop in order to change it for every found q(i). Before the for loop where i defined fuctions is working but when it goes in to for loop, program says that the fi1=f1(q(i)) fuctions is not defined. I tried many things but could not able to solve it. Also, I think my piecewise syntaxes are wrong. So if you can have a look at them it willl be geart! Thank you for your help !

clc;
fi1     = @(q)                        (atan(0.12*sin(q)/(0.4+0.12*cos(q))));
fi2     = @(q,fi1)                    (0.12*sin(q)/sin(fi1));
fi3     = @(fi1)                      (asin(-0.8*sin(fi1)/0.28));
fi4     = @(fi1,fi3)                  (0.92-0.8*cos(fi1)-0.28*cos(fi3));
xg3     = @(fi1)                      (0.4*cos(fi1));
yg3     = @(fi1)                      (0.4*sin(fi1));
g1      = @(q,fi1,fi2)                (0.12*cos(q-fi1)/fi2);
g2      = @(q,fi1)                    (0.12*cos(q+fi1));
g3      = @(g1,fi1,fi3)               (-0.8*cos(fi1)*g1/(0.28*cos(fi3)));
g4      = @(fi1,g1,fi3,g3)            (-0.8*sin(fi1)*g1-0.28*sin(fi3)*g3);
ug3     = @(fi1,g1)                   (-0.4*g1*sin(fi1));
vg3     = @(fi1,g1)                   (0.4*g1*cos(fi1));
g1dot   = @(q,g1,g2,fi1,fi2)          ((-2*g1*g2+0.12*sin(fi1-q))/fi2);
g2dot   = @(q,fi1,g1,fi2)             (-0.12*cos(q+fi1)+gi1^2*fi2);
g3dot   = @(g1,g1dot,g3,fi1,fi3)     
((0.8*(-sin(fi1)*g1^2+g1dot*cos(fi1))-0.28*sin(fi3)*g3^2)/(-0.28*cos(fi3)));
g4dot   =
@(g1,g1dot,g3,g3dot,fi1,fi3)(0.8*(cos(fi1)*g1^2+g1dot*sin(fi1))+0.28*(cos(fi3)*g3^2+g3dot*sin(fi3)));
ug3dot  = @(fi1,g1,g1dot)             (-0.4*g1dot*sin(fi1)-0.4*g1^2*cos(fi1));
vg3dot  = @(fi1,g1,g1dot)             (0.4*g1dot*cos(fi1)-0.4*g1^2*sin(fi1));
Fc      = @(q) piecewise((0.524<q)&&(q<2.618),500,(0.524>=q)|(q>=2.618),0);
T       = @(qdot)piecewise((0<=qdot/(2*pi/60))&&(qdot/(2*pi/60)<=1300),300*sin(-pi/650*(qdot/(2*pi/60))+650,qdot/(2*pi/60)>1300,-3,25*qdot/(2*pi/60)+4875));
J       = @(ug3,vg3,g4)                 (58.8+15*(ug3^2+vg3^2)+8*g4^2);
C       = @(ug3,ug3dot,vg3,vg3dot,g4,g4dot)(15*(ug3*ug3dot+vg3*vg3dot)+8*g4*g4dot);
Q       = @(vg3,g4,Fc,T)                (-15*9.81*vg3+Fc*g4+T);
dt      = 0.1;
t       = 0:dt:15;
t(1)    = 0;
q(1)    = 0;
qdot(1) = 0;
qddot   = @(q,qdot,t)((Q-C*qdot)/J);
for i=1:1:(length(t)-1)
    fi1= fil(q(i));
    fi2= fi2(q(i),fi1(i));
    fi3= fi3(fi1(i));
    fi4= fi4(fi1(i),fi3(i));
    xg3= xg3(fi1(i));
    yg3= yg3(fi1(i));
    g1 = g1(q(i),fi1(i),fi2(i));
    g2 = g2(q(i),fi1(i));
    g3 = g3(g1(i),fi1(i),fi3(i));
    g4 = g4(fi1(i),g1(i),fi3(i),g3(i));
    ug3= ug3(fi1(i),g(i));
    vg3= vg3(fi1(i),g1(i));
    g1dot= g1dot(q(i),g1(i),g2(i),fi1(i),fi2(i));
    g2dot= g2dot(q(i),fi1(i),g1(i),fi2(i));
    g3dot= g3dot(g1(i),g1dot(i),g3(i),fi1(i),fi3(i));
    g4dot= g4dot(g1(i),g1dot(i),g3(i),g3dot(i),fi1(i),fi3(i));
    ug3dot= ug3dot(fi1(i),g1(i),g1dot(i));
    vg3dot= vg3dot(fi1(i),g1(i),g1dot(i));
    Fc = Fc(q(i));
    T  = T(qdot(i));
    J  = J(ug3(i),vg3(i),g4(i));
    C  = C(ug3(i),ug3dot(i),vg3(i),vg3dot(i),g4(i),g4dot(i));
    Q  = Q(vg3(i),g4(i),Fc(i),T(i));
    qddot = qddot(qdot(i),J(i),C(i),Q(i));
    k11   = qdot(i);
    k12   = qddot(q(i),qdot(i),t(i));
    k21   = qdot(i)+dt/2*k12;
    k22   = qddot(q(i)+dt/2*k11,qdot(i)+dt/2*k12,t(i)+t/2);
    k31   = qdot(i)+dt/2*k22;
    k32   = qddot(q(i)+dt/2*k21,qdot(i)+dt/2*k22,t+dt/2);
    k41   = qdot(i)+dt*k32;
    k42   = qddot(q(i)+dt*k31,qdot(i)+dt*k32,t+dt);
    q(i+1)= q(i)+dt/6*(k11+2*k21+2*k31+k41);
    qdot(i+1)= qdot(i)+dt/6*(k12+2*k22+2*k32+k42);
end
figure(q-t);
plot(q,t);

[Shubham Gupta]

I have found some mistakes. Since, I can't access 'piecewise' I won't be able to help you further. Mistakes that I found are following :

1. You are getting confused between 'fil' and 'fi1' since l and 1 are hard to distinguish, make sure to double check your variable name.
2. g is not defined.
3. When you defined g2dot outside loop you made some mistake in variable name ( to be specific gi1 should be g1 imo )

[madhan ravi]

'piecewise' requires symbolic toolbox

Answer 1

[madhan ravi]

I am not going to check your code , download the file exchange below and adapt it to your needs:

https://www.mathworks.com/matlabcentral/fileexchange/29851-runge-kutta-4th-order-ode