I am trying to model a shock formation with two coupled ode's. I know the equations are right but the solver struggles with the shock and oscillates after it occurs. I have tried all of the solvers (ode45, ode23, ode23s etc.) to no avail. The best result I got was to use fixed time step 2nd order runge-kutta (downloaded from the internet ODE2) with a very small time step but this still gets oscillations. Can anybody recommend which options would help or which other solvers I can use?

Dear Dave!

> I am trying to model a shock formation with two coupled ode's. I know the equations are right but the solver struggles with the shock and oscillates after it occurs. I have tried all of the solvers (ode45, ode23, ode23s etc.) to no avail. The best result I got was to use fixed time step 2nd order runge-kutta (downloaded from the internet ODE2) with a very small time step but this still gets oscillations. Can anybody recommend which options would help or which other solvers I can use?

It is unlikely that a fixed step 2nd order Runge-Kutta produces "better" results than the adaptive 4/5-order Runge-Kutta of ODE45.
I know several system which tend to oscillate after a shock, e.g. the strings of a guitar, the air in a flute, the predator/prey relation of a Lotka-Volterra system, ...
Please explain what the "struggle" of the solver means. What does "the solver oscillate" mean? If this means an oscillation of the results, why are you convinced, that this is not the correct solution?

Jan

Dear Jan,

Thank you for your reply. I am trying to replicate a solution from a paper (so I know what it should look like). When I use anything (ode45, ode23 etc.) but the ODE1, ODE2, ... , ODE5. i.e. fixed time step solvers I get "Warning: Failure at t=2.335459e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
(7.105427e-15) at time t." So I have been trying to use the fixed time step solvers for very small values of t since all the other ode system solvers produce a similar error to the above. Are there any options that may help with the above error?

Thanks again

Is there any way I can post you a picture to show the comparison? It replicates the solution on either side of the shock but not the shock itself.

Dear Dave!

> Thank you for your reply. I am trying to replicate a solution from a paper (so I know what it should look like). When I use anything (ode45, ode23 etc.) but the ODE1, ODE2, ... , ODE5. i.e. fixed time step solvers I get "Warning: Failure at t=2.335459e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed
> (7.105427e-15) at time t." So I have been trying to use the fixed time step solvers for very small values of t since all the other ode system solvers produce a similar error to the above. Are there any options that may help with the above error?

What does this problem mean? The ODE solver detects a discontinuity at the time 2.335459. A fixed step solver runs through/over this discontiuity without detecting it - it works, but the is not a scientifically correct method.
What happens at this time point? Are you able to handle the discontinuity usign [Options.Events] ?

Good luck, Jan

Thank you again. I am trying the event option but the Matlab literature is confusing. Is there any chance you could give me some sort of a template. My code is quite simple

function dy = roberts(eta,y)

% define physical constants

omega=25;             % angular velocity at omega(r) (km /sec /kpc)
kappa=31.3;           % epicyclic frequency (km /sec /kpc)
omegapattern=13.5;    % angular pattern velocity (km /sec /kpc)
c=8.56;               % speed of sound (km /sec)
r=10;                 % radius of circular flow (kpc)
alpha=0.11667;        % angle of spiral inclination (radians)
A=72.92;              % amplitude (km /sec)^2

v0=r*(omega - omegapattern);
u0=alpha*r*(omega - omegapattern);

% function to be integrated

dy = zeros(2,1);

dy(1) = (y(1)/((y(1)^2) - (c^2)))*(2*omega*(y(2) - v0) + ((2*A)/(alpha*r))*sin(((2*eta)/(alpha*r))));
dy(2) = -((kappa^2)/(2*omega))*((y(1) - u0)/y(1));

I then run

etaspan = [0 3.41]

y0 = [13.41 115]

 [eta, y] = ode45(@roberts, etaspan, y0);

in the command window