Avoid divergent curves of ODE solutions

Hi all

I have solved a system of differential equations in an iterative way, as you can see from the graph I have several curves each of which is associated with a different combination of the iterative parameters. Again from the graph, it can be seen that some curves at a certain abscissa diverge; this is not acceptable so I would like to remedy it. Is it possible to create a new solution vector capable of avoiding the divergence as described below ?

I would like the solution not to diverge near the various cusp points and, on the contrary, to maintain the value immediately upstream of the latter, until the end of the domain. I'll explain: suppose we have the following solution values 1000 K,1020 K,2000 K,4000 K I would like instead: 1000 K,1020 K,1020 K,1020 K

Thank you for the help

Regards

T0=293;
delta_Thot=2500;%s 
%%%%%%%%%%%%%
Nc_Thot=2;
[ya,sa,lea]=size(Aexposehot);
for y=1:yV%passo,altezza tubo,Rextbobbin,lunghezza cavo
    for le=1:leV%lunghezza tubo equilibrio
                for sv=1:sV%spessore parete3
                  % for k=1:km%Mach 
                   for rf=1:length(Pressure_percentage_hot)
Rextbobb_cellhot(y)={Rextbobbin_fluxhot(y)};
Y0hot(y,le)={[T0,m0hot(y,le)]};
Dall_Th_hot ={[0:Nc_Thot:delta_Thot]};        
Lequilibrium_cellhot(le)={Lequilibriumhot(le)};        
Passo_cellhot(y)={Passo_fluxhot(y)}; 
hrad_cellhot(y)={Hradial_fluxhot(y)};
Spessore_parete3_cellhot(sv) ={Spessore_parete3hot(sv)};
Aexpose_cellhot(y,sv,le)={Aexposehot(y,sv,le)};
Volume_cell_Aisihot(le,y,sv)={Volume_Aisihot(le,y,sv)};
Volume_cell_gashot(y,le)={Volume_gashot(y,le)};
throat_cell_hot(y)={ThrtR_fluxhot(y)};
Pressure_end_chamber_hot_cell(rf)={Pressure_end_chamber_hot(rf)};
                   end
end    
end
  end


for iDom_Th = 1:numel(Dall_Th_hot)
   for leLE=1:numel(Lequilibrium_cellhot)
       for yY=1:numel(hrad_cellhot)
               for rF=1:numel(Pressure_end_chamber_hot_cell)
                  for sS=1:numel(Spessore_parete3_cellhot)

TSol_hot(iDom_Th,yY,leLE,sS,rF)={zeros(length(Dall_Th_hot{iDom_Th}),2)};

    end
    end 
    end
    end
   end 
opt_t_hot=odeset('NormControl','Refine','Stats','MaxStep');

for iDom_Th = 1:numel(Dall_Th_hot)
    tRange_hot = Dall_Th_hot{iDom_Th};
         for leLE=1:numel(Lequilibrium_cellhot)
                 for yY=1:numel(hrad_cellhot)
                       for rF=1:numel(Pressure_end_chamber_hot_cell)
                         for sS=1:numel(Spessore_parete3_cellhot)
 
                           Ltube_hot= Lequilibrium_cellhot{leLE};%L_cellhot{yY}; 
                           Leq_hot= Lequilibrium_cellhot{leLE} ;   
                           Vol_Aisi_hot=Volume_cell_Aisihot{leLE,yY,sS};
                           Aexp_hot=Aexpose_cellhot{yY,sS,leLE};
                           Hradial_hot=hrad_cellhot{yY};
                           ThroatradiuS_hot=throat_cell_hot{yY};
                           Pressure_hot_end_chamber=Pressure_end_chamber_hot_cell{rF};

                 if Vol_Aisi_hot>0&&Leq_hot>0 &&Aexp_hot>0 &&Ltube_hot>0 && ThroatradiuS_hot>0&&Pressure_hot_end_chamber>0  

    [tSol{iDom_Th,yY,leLE,sS,rF},TSol_hot{iDom_Th,yY,leLE,sS,rF}]=ode23t(@(t,Y) ...
            temperaturaheaterclassica(t,Y,Ltube_hot,Voltage1,Vol_Aisi_hot,Aexp_hot,ThroatradiuS_hot,Pressure_hot_end_chamber),tRange_hot,Y0hot{yY,leLE},opt_t_hot);
        if(any(TSol_hot{iDom_Th,yY,leLE,sS,rF}(:,1)>1600))
       %      disp([ ' condizione verificata'])
      %disp (iDom_Th)
      disp(yY)
      disp(leLE)
      %disp(sS)
      %disp(rF)
     break    
        end   
          
                         end
                      end
                    end
                 end
         end
end
    
%end
[tT,yT,leT,sT,rFT]=size(TSol_hot);

for t = 1:tT
    for y=1:yT
         for le=1:leT
                 for s=1:sT
                      for rf=1:rFT
     if TSol_hot{t,y,le,s,rf}(end,1)>0
        
                                
figure(1);set(gcf,'Visible', 'on')
plot(tSol{t,y,le,s,rf}(:,1),TSol_hot{t,y,le,s,rf}(:,1))
xlabel('time [s]')
ylabel('Teq-transition/turbulent [K]')    
hold on
end
end
end
end
end
end

[EldaEbrithil]

So in summary:

Is it possible to convert diverging curves to horizontal curves?

Why are these diverging curves present from a certain abscissa onwards?

[EldaEbrithil]

Any help?

Answer 1

[Julio Sanchez]

The divergence of an ode system could be due to two things: on one hand the system itself is devergent, this is due to the combination of parameter you are giving to the ode. On the other hand, it could be that the integrator is not suitable or correctly configurated to solve that syste. To test this you could use different integrators as: ode45, ode23, ode113, ode15s, ode23t, ode23b or ode23tb. I suggest ode15s is very robust.

[EldaEbrithil]

hi Julio thank you for the reply. I have just used ode15s and i see better results but some divergent curves is always present

[EldaEbrithil]

how can i filter those divergent curves?

[Julio Sanchez]

mmmm I don´t think those curves are divergent but thats the behavior of the system itself due to the combination of parameters. I suggest you to test individualy those problematic lines with other integrators, if the bechavior is the same then the system is divergent for that combination of parameters. To see which curves are divergent you have to test each one by one, so take small ranges for the parameters you defined and analyze the system.

[EldaEbrithil]

ok thank you very much Julio!!

[EldaEbrithil]

i have tried what you have said but i can not avoid to use thse parameters combination. As you can see, each curves reach a maximum then remain constant for a certain Delta t; after that:

1) remain constant till the end of the domain

2) decrease

3) increase

is there the possibility to cut those curves that exhibit 2) and 3) ? i don't want to write a code that allows me to avoid those curve but instead i need a code that preserves those curve till their "strange" behaviour

[Julio Sanchez]

Hi again! i want to ask you why do you want to cut those curves that exhibit a different behavior, is you already tried with different integrators and the answers is the same, then that is the inherent behavior of the system with those parameters. What I would do in your case is to separate the curves with different behavior as your already explained in the comment and analyze the parameter, you could find something interesting. If you just want to "eliminate" the "strange" behavior then you could just put an if condition so that the returned value of the ode is zero, then the system will remain in that state.

[EldaEbrithil]

Hi Julio

my problem is that i am not able to set up an if condition suitable for eliminating those curves:

the condition

if(any(TSol_hot{iDom_Th,yY,leLE,sS,rF}(:,1)>1600))

is not useful because as you can see from the plots it doesn't cut the curve that exceed T>1600

moreover i obtain a lot of these warnings:

Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.637979e-12) at time t. 
> In ode15s (line 668)

for this reason i started my topic thinking that those curves may be related with unacceptable behaviours of the system

[EldaEbrithil]

More in detail:

i need only those solutions which exibhit an horizontal terminal part

[Julio Sanchez]

Hi! then, the only easy way (without to much code) to see which lines have the undesired behavior, is again to test the combinations of parameters one by one and discard the divergent ones.

[EldaEbrithil]

Hi Julio

is there the possibility to go though each step inside a for loop? In this way i could be able to test each curve one by one

[Julio Sanchez]

That's a good solution! you can put integrator inside the for loop like this:

storage = {}
for i=1:length of parameter
    storage{i} = integrator(parameter(i))
end

Then, as the solutions are inside the storage cell, you can choose wich ones reach steady state and only plot those.

[EldaEbrithil]

something like that?

storage = {};
for t = 1:tT
    for y=1:yT
         for le=1:leT
           for s=1:sT
               for rf=1:rFT
                    for pE=1:pET
    storage{t,y,le,s,rf,pE} =  TSol_hot(t,y,le,s,rf,pE);
                          end
                      end
           end
         end
    end
end
 storage2=cell2mat(storage);

however TSol_hot is constituted by 1x1 cells each one is {2251×2 double}. It gives me error on storage2

[Julio Sanchez]

The problem here is that you defined the sorage cell as a multidimentional space with t, y, le,rf and pE as positions inside the cell. You could just do this for your problem:

storage(i) = {[t,y,le,s,rf,pE],Tsol_hot(t,y,le,s,rf,pE)}

That way each position in the storage cell has two components: a vector of parameters and a matrix with the solution of the ode system. The counter i length must be then the sum of the lengths of tT, yT, leT, sT, rfT and peT.