Ode45 not able to plot anything

1 view (last 30 days)
EldaEbrithil
EldaEbrithil on 18 Jun 2021
Commented: EldaEbrithil on 18 Jun 2021
Hi all
i have some trouble with this iterative differential problem:
Teta0=0;
dTeta0=0;
Nc=0.1;
tempo=10;
for k=1:length(dTeta0)
for i=1:length(Lung_finale)
Y0(k)={[dTeta0(k),Teta0]};%Initial condition vector
Time={[1:Nc:tempo]};
kd_cell(i)={kd_finale(i)};
Ix_cell(i)={Ix_finale(i)};
D_cell(i)={D_finale(i)};
h_cell(i)={H_finale(i)};
Ixd_cell(i)={Ix_demihull_finale(i)};
Awd_cell(i)={Awd_finale(i)};
L_cell(i)={Lung_finale(i)};
Larg_cell(i)={Larg_finale(i)};
Cw_cell(i)={Cw_finale(i)};
T_cb_cell(i)={T_cb_finale(i)};
Awc_cell(i)={Awc_finale(i)};
V_cell(i)={V_finale(i)};
B_cell(i)={B_finale(i)};
CB_cell(i)={CB_finale(i)};
end
end
for iDom = 1:numel(Time)
for iInitial=1:numel(Y0)
for i=1:numel(L_cell)
YSol(iDom,iInitial,i)={zeros(length(Time{iDom}),2)};
end
end
end
for iDom = 1:numel(Time)
tRange = Time{iDom};
for iInitial=1:numel(Y0)
for i=1:numel(L_cell)
kd_ode=kd_cell{i};
Ix_ode=Ix_cell{i};
D_ode=D_cell{i};
h_ode=h_cell{i};
Ixd_ode=Ixd_cell{i};
Awd_ode=Awd_cell{i};
L_ode=L_cell{i};
Larg_ode=Larg_cell{i};
Cw_ode=Cw_cell{i};
T_ode=T_cb_cell{i};
Awc_ode=Awc_cell{i};
V_ode=V_cell{i};
B_ode=B_cell{i};
CB_ode=CB_cell{i};
[tSol{iDom,iInitial,i},YSol{iDom,...
iInitial,i}]=ode45(@(t,Y)Rollio_ED_catamarano(t,Y,kd_ode,Ix_ode,D_ode,h_ode,Ixd_ode,Awd_ode,gamma,L_ode,Larg_ode,Cw_ode,CB_ode,rho,...
T_ode,Awc_ode,B_ode,V_ode),tRange,Y0{iInitial});
end
end
end
for iDom = 1:numel(Time)
tRange = Time{iDom};
for iInitial=1:numel(Y0)
for i=1:numel(L_cell)
figure(1);set(gcf,'Visible', 'on')
if isreal(YSol{iDom,iInitial,i}(:,2))
plot(tSol{iDom,iInitial,i}(:,1),YSol{iDom,iInitial,i}(:,2))
xlabel('Time [s]')
ylabel('Roll angle [rad]')
hold on
end
end
end
end
As u can see, i want to solve iteratively an ODE problem in a fixed time domain.
The differential equation is over there:
function dYdt=Rollio_ED_catamarano(t,Y,kd_ode,Ix_ode,D_ode,...
h_ode,Ixd_ode,Awd_ode,gamma,L_ode,Larg_ode,Cw_ode,CB_ode,rho,...
T_ode,Awc_ode,B_ode,V_ode)
m=2.5;
pc=(Awc_ode*kd_ode^2)/(V_ode*h_ode);
v=0.1164;
Omega_ode=(0.5)*t;%<------ time dependent
lambda_w_ode=1.56*(2*pi/Omega_ode)^2;%<------ time dependent
k=(2*pi/lambda_w_ode);
k_bar=k/Larg_ode;
n=1+2.5*k_bar^m;
m_teta=1-0.3*k_bar^2;
X=CB_ode/Cw_ode;
X_tetat_0=1-((6*X^3)/(1+X)*(1+2*X))*k*T_ode+(((1.5*X^3)/(2-X)*(2+X))*(k*T_ode)^2)-((X^3)/(3*(3-2*X)*(3-X)))*(k*T_ode)^3;%Yung pag 232
X_tetaB=m_teta*(1-sqrt(Cw_ode)*(B_ode/lambda_w_ode)^2);
Xc=1+(0.5/(0.5+0.75+sqrt(k_bar)));
X_tetat=Xc*X_tetat_0;
X_teta=X_tetaB*X_tetat;
X_ziT_0=1+X*k*T_ode+(X/(2*(2-X)))*((k*T_ode)^2)-(X/(6*(3-2*X)))*(k*T_ode)^3;
X_zib=1-1.73*Cw_ode*((Larg_ode*(1+(n/(n-k_bar))))/(lambda_w_ode))^2;
X_zi=X_zib*Xc*X_ziT_0;
f=X_teta/X_zi;
Lamb_33d=(0.85*pi/4)*(gamma/9.81)*L_ode*(Larg_ode^2)*((Cw_ode^2)/(1+Cw_ode));
Delta_lamb_33=0.21*rho*pi*T_ode*L_ode*((Larg_ode^3)*Cw_ode^3)/((kd_ode^2)*(1+Cw_ode)*(1+2*Cw_ode));
X1_curva_ode = (4.28e-03)*(L_ode/lambda_w_ode)^4 - (7.01e-02)*(L_ode/lambda_w_ode)^3 +...
(4.15e-01)*(L_ode/lambda_w_ode)^2 - 1.07*(L_ode/lambda_w_ode) + 1.11;%<------ time dependent
X2_curva_ode=+(3.79e+01)*((T_ode/lambda_w_ode)^4)- 5.85e+01*((T_ode/lambda_w_ode)^3) +...
3.34e+01*((T_ode/lambda_w_ode)^2) - 8.83e+00*((T_ode/lambda_w_ode)) + 9.94e-01 ;%<------ time dependent
Lamb_33=2*Lamb_33d+2*Delta_lamb_33;
Lambda_zi=1.6*Lamb_33d;
Lamb_44=2.5*Lambda_zi*kd_ode^2;
eps=(2*Lambda_zi)/((D_ode/9.81)+2*Lambda_zi);
Omega_h=sqrt((2*gamma*Awc_ode)/((D_ode/9.81)+2*Lambda_zi));
Omega_r=sqrt(D_ode*h_ode/(Ix_ode+Lamb_44));
v_primo=(v*kd_ode^2)/((Ix_ode+Lamb_44)*Omega_r);
v_zi=0.5*rho*((Omega_ode^3)/9.81)*(Awd_ode^2)*(X2_curva_ode)*X1_curva_ode;
k_primo=Larg_ode/4;
Lamb_44d=2*Lamb_33d*k_primo^2;
v_teta=0.1*sqrt((D_ode*h_ode/2)*(Ixd_ode+Lamb_44d));
N_teta=v_teta+v_zi*kd_ode^2;
vc=N_teta/(Ix_ode+Lamb_44);
hw=0.17*lambda_w_ode^(3/4);%<------ time dependent
r=0.5*hw;
d=(cos(k*kd_ode))/(sin(k*kd_ode));
Lambda_teta=(Lamb_44-2*Lambda_zi*kd_ode^2)/2;
q_lambda=1+(Lambda_teta/(Lambda_zi*kd_ode^2));
q_v=1+(v_teta/(v_zi*kd_ode^2));
alfa=1;
Teta_dot=Y(1);
Teta=Y(2);
dTeta_dotdt=-2*vc*Teta_dot-(Omega_r^2)*Teta+X_zi*k*r*(Omega_r^2)*(sin(k*kd_ode)/k*kd_ode)*...
(sqrt((1+d^2)*(f*k*kd_ode-q_lambda*pc*eps*(Omega_ode^2)/(Omega_h^2))+4*q_v*(v_primo^2)*((Omega_ode^2)/(Omega_r^2))))*sin(Omega_ode*t+alfa);
dYdt=[dTeta_dotdt;Teta_dot];
end
I have noticed that if i decrease Omega (for example setting it to (1/100)*t) i was able to plot something, however these Omega values are too small for my analysis. Decresing Omega increases hw and lambda_w_ode. With bigger Omega hw and lambda_w_ode decrease, may it be a problem in the solving process? Do you have any idea behind this fact?
Thank you very much for the help
Regards!!
  4 Comments
EldaEbrithil
EldaEbrithil on 18 Jun 2021
Tank you very much for the hints Jan, they have been useful!!

Sign in to comment.

Answers (0)

Categories

Find more on Parallel for-Loops (parfor) in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by