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

Question

DJ V 2024-9-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152635-matlab-ode-4-5-runge-kutta-problem-solver-can-t-get-it-to-work-for-one-problem-set

编辑： Sam Chak 2024-9-13

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

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

请先登录，再进行评论。

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

Answer 1

Sam Chak 2024-9-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2152635-matlab-ode-4-5-runge-kutta-problem-solver-can-t-get-it-to-work-for-one-problem-set#answer_1516335

在 MATLAB Online 中打开

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

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

请先登录，再进行评论。

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

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

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

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

回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

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