Results 1 to 3 of 3

Math Help - Matlab ODE solver issues

  1. #1
    Junior Member
    Joined
    Oct 2009
    Posts
    31

    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?
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by davefulton View Post
    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
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Junior Member
    Joined
    Oct 2009
    Posts
    31
    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
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Complex Solver
    Posted in the Algebra Forum
    Replies: 4
    Last Post: April 10th 2011, 02:24 AM
  2. Numeric Solver in Matlab
    Posted in the Math Software Forum
    Replies: 1
    Last Post: April 25th 2010, 07:30 PM
  3. Cubic Solver
    Posted in the Math Software Forum
    Replies: 0
    Last Post: February 10th 2009, 03:10 PM
  4. Anybody interested in a Quadratic Solver?
    Posted in the Math Software Forum
    Replies: 52
    Last Post: January 11th 2008, 05:49 AM
  5. Differential Equation Solver.
    Posted in the Calculators Forum
    Replies: 0
    Last Post: March 28th 2005, 04:41 PM

Search Tags


/mathhelpforum @mathhelpforum