Reduce computational time for ODE solver

Hi all

i have this code

for iDom_hot = 1:numel(Dall_hot)
    xRange_hot = Dall_hot{iDom_hot};
for iInitial_hot = 1:numel(Y0_Lhot)

            Mach_hot=Mach_cell_hot{iInitial_hot};
            Tw_hot=Tnew_cell_hot{iDom_hot};
            Ch_hot=ChL_cell_hot{iDom_hot,iInitial_hot};
            Fc_hot=fL_cell_hot{iDom_hot,iInitial_hot};   
            Hrrj_hot=hrrj_cell_hot{iDom_hot};
            Aeff_hot=Aeff_cell_hot{iDom_hot,iInitial_hot};
            Perim_hot=Perim_cell_hot{iDom_hot,iInitial_hot};
            L_Ength_hot=L_cell_hot{iDom_hot};
            Pend_chamb_hot=Pressure_ED_chamber_cell_hot{iDom_hot};
         if Ch_hot>0&&Fc_hot>0&&Hrrj_hot>0&&Tw_hot>0 &&Mach_hot>0&&Perim_hot>0&&Aeff_hot>0&&L_Ength_hot>0&&Pend_chamb_hot>0
             
    [xSol_hot{iDom_hot,iInitial_hot},YSol_L_hot{iDom_hot,iInitial_hot}]=ode23(@(x,Y)...
            chambesinglebobbTWVARIABILEHOT(x,Y,Ch_hot,Aeff_hot,Perim_hot,Fc_hot,Tw_hot),xRange_hot,Y0_Lhot{iInitial_hot});  
        if(any(YSol_L_hot{iDom_hot,iInitial_hot}(:,3)>0.8))
      disp([ ' condizione verificata'])
      disp (iDom_hot)
      disp(iInitial_hot)
      
     break    
   end     
         end
         
   
   
            
end
           
end

that solve iteratively an ODE system:

function dYdx=tubosinglebobb(x,Y,Ch_hot,Aeff_hot,Perim_hot,Fc_hot,Tw_hot)
gamma=1.667;
Dexttantalum=0.001500000000000;
Tt=Y(1);
Pt=Y(2);
M=Y(3);

dTtdx=Ch_hot*(Tw_hot-Tt)*(Perim_hot/Aeff_hot);
dPtdx=Pt*((-gamma*((M^2)/2)*Fc_hot*(Perim_hot/Aeff_hot))-...
    (gamma*((M^2)/2)*dTtdx*(1/Tt)));
dMdx=M*(((1+((gamma-1)/2)*M^2)/(1-M^2))*((gamma*(M^2)*Fc_hot*...
    Perim_hot)/(2*Aeff_hot))+...
    ((0.5*((gamma*M^2)+1))*(1+((gamma-1)/2)*M^2)/(1-M^2))*...
    (Ch_hot*(Tw_hot-Tt)*Perim_hot)/(Aeff_hot*Tt));

dYdx=[dTtdx;dPtdx;dMdx];
end

the computational time for solving the ODE system is very very large (days) due to ( i think) one problem. In particular:

numel(Dall_ho) is more or less equal to 700, while numel(Y0_Lhot) is equal to 30, those domains dimensions should not be problematic, the real problem is on the oscillation of the solution YSol_L_hot when its third column reach value around 1. I have prepared a break for this eventuality but doesn't seem to reduce computational time!! I don't really know which is the problem, could it be a syntax error?

Thank you all for the help

Regards

Answer 1

[Bjorn Gustavsson]

Have you tried the other ODE-integrating functions? From your problem description it sounds like your ODEs are stiff - for that ode23 is not the optimal function, perhaps you get much better success with ode15s or oed23s that both are designed for stiff ODEs. One further note is that to me it seems like the temperature-derivative does not depend on the other 2 parameters, perhaps you can refactor the problem to first integrate that ODE and then use that solution when integrating the other 2.

[EldaEbrithil]

Thank you for the reply

i have tried ode23s but without success, the computational time it is the same. I also have used ode15s, this solver is extremly rapid but it gives me this warning numerous time:

Warning: Failure at t=1.133336e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

> In ode15s (line 668)

[Bjorn Gustavsson]

Well if none of the odeNN-functions works, perhaps you can figure out what's happening when ode15s sends the warning if you set the debug-flag to stop when warning:

dbstop if warning

or

dbstop in ode15s if t>1.13333e0

Most of the time I run into this it is because of a simplification/approximation or assumption about the solution in my ode-function is breaking down in some way or the other.

[EldaEbrithil]

Thank you very much!!

[EldaEbrithil]

paradoxically ode45 solve the problem without warnings and in short times, but in a not accurate way (for example the Mach number M becomes negative...this is impossible from the physics of equations), i really don't know why this fact happens...

[Bjorn Gustavsson]

If you're dealing with some kind of transition (bounces, phase-changes or shocks going between sub and supersonic) you might want to capture that with some extra care, best way of doing that is to use the event-handling that the odeNN-functions have, check the ballode for how to use the 'Events' handling.

[EldaEbrithil]

i have inserted thsi line using ode45 as sovler:

opt=odeset('RelTol',1e-5,'AbsTol',1e-7);

this line should increase accuracy for ode 45 but it only increases the time....

[EldaEbrithil]

Moreover ode 45 is very less accurate than ode 23, the results are inaccetable... i seriuously don't know how to do for speeding up the process keeping a moderate accuracy

[Bjorn Gustavsson]

Since you're modelling some physical system, where you have some assumptions and some physical constraints on what makes a valid solution. You have to make some analysis of what happens with your solutions (M < 0) and also make some checks on what happens around those periods. For handling things like this you can use the options to force (some components of) the solution to be positive: ode-options, and to handle shocks and transitions you can use the even-handling: event-location.