Matlab ODE solver issues

Oct 2009
31
0
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?
 

CaptainBlack

MHF Hall of Fame
Nov 2005
14,975
5,273
erehwon
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
 
Oct 2009
31
0
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