Ode45 not able to plot anything

2 Ansichten (letzte 30 Tage)
EldaEbrithil
EldaEbrithil am 18 Jun. 2021
Kommentiert: EldaEbrithil am 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 Kommentare
EldaEbrithil
EldaEbrithil am 18 Jun. 2021
Bearbeitet: EldaEbrithil am 18 Jun. 2021
I am sorry, you are right...
the problem is that I am currently not able to obtain a plot for my differential problem, I suppose this is due the magnitude of some time dependent parameter inside my equation. What I have found is that after I had changed my time interaval a little (Time 0.01:0.1:1) I was able to obtain some plot:
Changing further the initial/last Time extremes leads to no-graph. It is interesting to underline that, when I had set Omega_ode to very small values (for example using the equation Omega_ode=0.001*t) I was able to plot my graph, however these Omega_ode values are not the ones I need for my analysis.
I'm trying to figure out why this happens
EldaEbrithil
EldaEbrithil am 18 Jun. 2021
Tank you very much for the hints Jan, they have been useful!!

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by