# Thread: Matlab ODE solver issues

1. ## Matlab ODE solver issues

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?

2. Originally Posted by davefulton
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?
Could you post the equations?

(the sharpness of the shock can be a problem. All the numerical methods effectively assume smooth solutions on some scale, so you may be seeing an analogue of Gibbs phenomenon, but then I'm not an expert on modelling shocks)

CB

3. HI Captain Black,

Thank you for your reply. I missed out a periodicity condition so the equations can be solved with the bvp4c solver. However I get the same problem. The Matlab code is

Code:
function BVPattempt

r=10;                 % radius of circular flow (kpc)
alpha=0.11667;        % angle of spiral inclination (radians)

solinit = bvpinit(linspace(0,alpha*r*pi,1000),[13.4171 115]); % not sure about the last argument, its a constant guess for the solution across the grid

sol = bvp5c(@f,@ex1bc,solinit);

% The solution at the mesh points
x = sol.x;
y = sol.y;

figure;
plot(x,y(1,:)')
hold on
plot(x,y(2,:)')
legend('u','v');
xlabel('eta');
ylabel('u/v');
title('None');
%---------------------------------------------------------------

function dudy = f(y,u)
dudy = zeros(2,1);

% 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);

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

%---------------------------------------------------------------

% Not sure how to impose two boundary conditions

function res = ex1bc(ya,yb)

%res = [ ya(1) - 13.4171
%yb(1) - 13.4171];
res = [ ya(1) - 13.4171
yb(1) - 13.4171];
Thanks again