Matlab ODE 4.5 Runge Kutta problem solver. Can't get it to work for one problem set.

I can't get the problem solver to work for myFun4. I don't understand why. It appears like the prior three myFun1, myFun2, and myFun3. The file is below.

type RG1Combined.m 

theta = 0;
phi = 0;
psi = 0;
Ctheta=cos(theta);
Cphi = cos(phi);
Cpsi = cos(psi);
Stheta = sin(theta);
Sphi = sin(phi);
Spsi = sin(psi);
Jx = 0.8244;
Jy= 1.135;
Jz = 1.759;
Jxz = 0.1204;
G = Jx*Jz-Jxz^2;
G1 = (Jxz*(Jx+Jy+Jz)/G);
G2 = (Jz*(Jz-Jy)+Jxz^2)/G;
G3 = Jz/G;
G4 = Jxz/G;
G5 = (Jz - Jx)/Jy;
G6 = Jxz/Jy;
G7 = ((Jx - Jy)*Jx +Jxz^2)/G;
G8 = Jx/G;
p = 0;
q = 0;
r = 0;
u = 0;
v = 0;
w = 0;
%m = mass.
m = 11;

tRange = [0, 10];
[tSol,uvwAndrSol] = ode45(@(t,uvwAndr)myFun1(t,uvwAndr,Ctheta,Cphi,Cpsi,Stheta,Sphi,Spsi),tRange,[1,1,1]);
length(tSol)
tSol
uvwAndrSol
uvwAndrSol(20,3)


[tSol,xyzAndrSol] = ode45(@(t,xyzAndr)myFun2(t,xyzAndr,p,q,r,u,v,w,m), tRange, [1,1,1]);
length(tSol)
tSol
xyzAndrSol
xyzAndrSol(20,3)

[tSol,pqrAndrSol] = ode45(@(t,pqrAndr)myFun3(t,pqrAndr,p,q,r,phi,psi,theta), tRange, [1,1,1]);
length(tSol)
tSol
pqrAndrSol
pqrAndrSol(20,3)

[tSol,lmnAndrSol] = ode45(@(t,lmnAndr)myFun4(t,lmnAndr,p,q,r,l,m,n,Jx,G1,G2,G3,G4,G5,G6,G7,G8), tRange, [1,1,1]);
length(tSol)
tSol
lmnAndrSol
lmnAndrSol(20,3)

function d_uvwAndr_dt= myFun1(t,uvwAndr,Ctheta,Cphi,Cpsi,Stheta,Sphi,Spsi)
%
d_uvwAndr_dt = zeros(3,1);
d_uvwAndr_dt(1) = Ctheta*Cpsi*Sphi*Stheta*Spsi*uvwAndr(1)-Cphi*Spsi*Cphi*Stheta*Cpsi*uvwAndr(2)+Sphi*Spsi*uvwAndr(3);
d_uvwAndr_dt(2) = Ctheta*Spsi*Stheta*Sphi*Spsi*uvwAndr(1)+Cphi*Cpsi*Cphi*Stheta*Spsi*uvwAndr(2)-Sphi*Cpsi*uvwAndr(3);
d_uvwAndr_dt(3) = -Stheta*uvwAndr(1)+Sphi*Ctheta*uvwAndr(2)+Cphi*Ctheta*uvwAndr(3);
%% 
end


function d_xyzAndr_dt=myFun2(t,xyzAndr, p, q, r, u,v,w,m)
%
d_xyzAndr_dt = zeros(3,1);
d_xyzAndr_dt(1) = (r*v-q*w) + (5*t)/m;
d_xyzAndr_dt(2) = (p*w-r*u) + (t)/m;
d_xyzAndr_dt(3) = (q*u-p*v) + (3*t)/m;
%
end

function d_pqrAndr_dt=myFun3(t,pqrAndr,p,q,r,phi,psi,theta)
d_pqrAndr_dt = zeros (3,1);
d_pqrAndr_dt(1) = 1*p + sin(psi)*tan(theta)*q + cos(psi)*tan(theta)*r;
d_pqrAndr_dt(2) = 0*p + cos(phi)*q - sin(phi)*r;
d_pqrAndr_dt(3) = 0*p + sin(phi) / cos(theta) * q +cos(phi) / cos(theta) * r;
end


function d_lmnAndr_dt=myFun4(t,lmnAndr,p,q,r,l,m,n,Jx,G1,G2,G3,G4,G5,G6,G7,G8)
d_lmnAndr_dt = zeros (3,1);
d_lmnAndr_dt(1) = G1*p*q-G2*q*r+G3*l+0/Jx+G4*n;
d_lmnAndr_dt(2) = G5*p*r-G6*(p^2-r^2)+0*l+m/Jx+0*n;
d_lmnAndr_dt(3) = G7*p*q -G1*q*r+G4*l+0/Jx+G8*n;
end

Answer 1

[Sam Chak]

Hi @DJ V

You forgot to supply the values for l and n.

theta = 0;
phi = 0;
psi = 0;
Ctheta=cos(theta);
Cphi = cos(phi);
Cpsi = cos(psi);
Stheta = sin(theta);
Sphi = sin(phi);
Spsi = sin(psi);
Jx = 0.8244;
Jy= 1.135;
Jz = 1.759;
Jxz = 0.1204;
G = Jx*Jz-Jxz^2;
G1 = (Jxz*(Jx+Jy+Jz)/G);
G2 = (Jz*(Jz-Jy)+Jxz^2)/G;
G3 = Jz/G;
G4 = Jxz/G;
G5 = (Jz - Jx)/Jy;
G6 = Jxz/Jy;
G7 = ((Jx - Jy)*Jx +Jxz^2)/G;
G8 = Jx/G;
p = 0;
q = 0;
r = 0;
u = 0;
v = 0;
w = 0;
%m = mass.
m = 11;

%% User-supplied parameters
l   = 1;
n   = 1;

tRange = [0, 10];
[tSol,uvwAndrSol] = ode45(@(t,uvwAndr)myFun1(t,uvwAndr,Ctheta,Cphi,Cpsi,Stheta,Sphi,Spsi),tRange,[1,1,1]);
length(tSol)

ans = 45

tSol

tSol = 45×1

0 0.0502 0.1005 0.1507 0.2010 0.4218 0.6426 0.8635 1.0843 1.3343

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

uvwAndrSol

uvwAndrSol = 45×3

1.0000 1.0000 1.0000 1.0000 1.0000 1.0515 1.0000 1.0000 1.1057 1.0000 1.0000 1.1627 1.0000 1.0000 1.2226 1.0000 1.0000 1.5248 1.0000 1.0000 1.9016 1.0000 1.0000 2.3714 1.0000 1.0000 2.9575 1.0000 1.0000 3.7980

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

uvwAndrSol(20,3)

ans = 46.2637

[tSol,xyzAndrSol] = ode45(@(t,xyzAndr)myFun2(t,xyzAndr,p,q,r,u,v,w,m), tRange, [1,1,1]);
length(tSol)

ans = 41

tSol

tSol = 41×1

0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xyzAndrSol

xyzAndrSol = 41×3

1.0000 1.0000 1.0000 1.0142 1.0028 1.0085 1.0568 1.0114 1.0341 1.1278 1.0256 1.0767 1.2273 1.0455 1.1364 1.3551 1.0710 1.2131 1.5114 1.1023 1.3068 1.6960 1.1392 1.4176 1.9091 1.1818 1.5455 2.1506 1.2301 1.6903

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xyzAndrSol(20,3)

ans = 4.0767

[tSol,pqrAndrSol] = ode45(@(t,pqrAndr)myFun3(t,pqrAndr,p,q,r,phi,psi,theta), tRange, [1,1,1]);
length(tSol)

ans = 41

tSol

tSol = 41×1

0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

pqrAndrSol

pqrAndrSol = 41×3

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

pqrAndrSol(20,3)

ans = 1

[tSol,lmnAndrSol] = ode45(@(t,lmnAndr)myFun4(t,lmnAndr,p,q,r,l,m,n,Jx,G1,G2,G3,G4,G5,G6,G7,G8), tRange, [1,1,1]);
length(tSol)

ans = 53

tSol

tSol = 53×1

0 0.0038 0.0075 0.0113 0.0151 0.0339 0.0527 0.0715 0.0904 0.1845

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

lmnAndrSol

lmnAndrSol = 53×3

1.0000 1.0000 1.0000 1.0049 1.0502 1.0025 1.0099 1.1005 1.0050 1.0148 1.1507 1.0074 1.0197 1.2010 1.0099 1.0444 1.4521 1.0223 1.0690 1.7033 1.0347 1.0936 1.9545 1.0471 1.1183 2.2057 1.0595 1.2415 3.4616 1.1214

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

lmnAndrSol(20,3)

ans = 2.4589

function d_uvwAndr_dt= myFun1(t,uvwAndr,Ctheta,Cphi,Cpsi,Stheta,Sphi,Spsi)
    %
    d_uvwAndr_dt = zeros(3,1);
    d_uvwAndr_dt(1) = Ctheta*Cpsi*Sphi*Stheta*Spsi*uvwAndr(1)-Cphi*Spsi*Cphi*Stheta*Cpsi*uvwAndr(2)+Sphi*Spsi*uvwAndr(3);
    d_uvwAndr_dt(2) = Ctheta*Spsi*Stheta*Sphi*Spsi*uvwAndr(1)+Cphi*Cpsi*Cphi*Stheta*Spsi*uvwAndr(2)-Sphi*Cpsi*uvwAndr(3);
    d_uvwAndr_dt(3) = -Stheta*uvwAndr(1)+Sphi*Ctheta*uvwAndr(2)+Cphi*Ctheta*uvwAndr(3);
    %% 
end

function d_xyzAndr_dt=myFun2(t,xyzAndr, p, q, r, u,v,w,m)
    %
    d_xyzAndr_dt = zeros(3,1);
    d_xyzAndr_dt(1) = (r*v-q*w) + (5*t)/m;
    d_xyzAndr_dt(2) = (p*w-r*u) + (t)/m;
    d_xyzAndr_dt(3) = (q*u-p*v) + (3*t)/m;
    %
end

function d_pqrAndr_dt=myFun3(t,pqrAndr,p,q,r,phi,psi,theta)
    d_pqrAndr_dt = zeros (3,1);
    d_pqrAndr_dt(1) = 1*p + sin(psi)*tan(theta)*q + cos(psi)*tan(theta)*r;
    d_pqrAndr_dt(2) = 0*p + cos(phi)*q - sin(phi)*r;
    d_pqrAndr_dt(3) = 0*p + sin(phi) / cos(theta) * q +cos(phi) / cos(theta) * r;
end

function d_lmnAndr_dt=myFun4(t,lmnAndr,p,q,r,l,m,n,Jx,G1,G2,G3,G4,G5,G6,G7,G8)
    d_lmnAndr_dt = zeros (3,1);
    d_lmnAndr_dt(1) = G1*p*q-G2*q*r+G3*l+0/Jx+G4*n;
    d_lmnAndr_dt(2) = G5*p*r-G6*(p^2-r^2)+0*l+m/Jx+0*n;
    d_lmnAndr_dt(3) = G7*p*q -G1*q*r+G4*l+0/Jx+G8*n;
end