Ode45 not able to plot anything

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!!

[EldaEbrithil]

Any idea?

[Jan]

Please bump your question after at least 24 hours only. Thanks.

Instead of posting a meaningless comment, care for explaining, what the problem is. Currently all we know is, that you have "some troubles". For posting a solution, it is required to know, which troubles you have. preferrably as exact as possible.

Some hints: This is just a waste of time:

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

This is no pre-allocation, but you let the array YSol grow iteratively. Defining a cell on the right side and copying it to the left side is a waste of time also. Better:

YSol = cell(numel(Time), numel(Y0), numel(L_cell));  % Pre-allocation
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);
          % The curly braces on the left side!
      end
   end  
end

But the contents of the cell is overwritten, so a pre-allocation is not useful at all. The best method is a simple:

YSol = cell(numel(Time), numel(Y0), numel(L_cell));  % Pre-allocation

This is not efficient also:

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)};
      ...          

Again, use the curly braces on the left side, not on the right. But this is worse: You redefine Y0{k} length(Lung_finale) times, and kd_cell{i} length(dTeta0) times, and Time even length(Lung_finale)*length(dTeta0) times with the same value. Avoid repeated work to reduce the run time and the confusion level.

[EldaEbrithil]

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]

Tank you very much for the hints Jan, they have been useful!!