How to use ode45 (or other) to solve equations without taking too big of a step along interval

I am trying to use ode23s to solve a ski jumper dynamics problem. I want to calculate the speed along path with given initial conditions. First, I create x and y variables to define the position of an approach, transition, and takeoff. Then I use ode45 to solve for the velocity. The output data skips over important coordinates, namely from x=horizantal position=53 to x=59, which is the entire transition and where the velocity change is most dynamic and where I most want data. How do I avoid this?

MATLAB Code

% File calculates the rider position along the x axis the velocity along
% the slope(i.e. sdot) and the acceleration along the slope(sdoubledot) for
% any point within a time interval tspan.
%
% Ap = approach
% Tr = transition
% To = takeoff
% L1 = inrun slope length
% L2 = Takoff slope length
% R = radius of transition

format compact
global xInrun yInrun dydxInrun kurvature
global m CdA mu eta g

L1 = 60;
L2 = 3;
R = 9;
thetaAp = 24*pi/180;
thetaTO = 15*pi/180;
xTrStart = L1*cos(thetaAp);
yTrStart = -L1*sin(thetaAp);
xToStart = xTrStart + R*sin(thetaAp) + R*sin(thetaTO);
yToStart = yTrStart - (R*cos(thetaTO) - R*cos(thetaAp));
xTOEnd = xToStart + L2*cos(thetaTO);
yTOEnd = yToStart + L2*sin(thetaTO);

% (x,y) data for inrun
numpnts = 100;
xAp = linspace(0,xTrStart,numpnts);
yAp = xAp.*-tan(thetaAp);

% (x,y) data for transition with radius R centered at (xTrOrigin,
% yTrOrigin)
xTrOrigin = xTrStart + R*sin(thetaAp);
yTrOrigin = yTrStart + R*cos(thetaAp);
xTr = linspace(xTrStart,xToStart,numpnts);
yTr = yTrOrigin - sqrt(R*R-(xTr-xTrOrigin).*(xTr-xTrOrigin));

% (x,y) data for takeoff
xTo = linspace(xToStart,xTOEnd,numpnts);
yTo = yToStart + (xTo-xToStart).*tan(thetaTO);

% combining data into one vector for velocity analysis
xInrun_raw = [xAp(1:99), xTr(1:99), xTo];
yInrun_raw = [yAp(1:99), yTr(1:99), yTo];
xInrun = linspace(0,xInrun_raw(end), 600);
yInrun = interp1(xInrun_raw,yInrun_raw,xInrun);

% calculate inrun speed from start to takeoff

ft2m=12/39.37;      % feet to meter conversion factor
m2ft=1/ft2m;        % meters to feet conversion factor
deg2rad=pi/180;   % degree to radian conversion factor
mps2mph=2.237;      % meters per second to miles per hour conversion
mph2mps=1/mps2mph;  % miles per hour to meters per second conversion

figure(1)
plot(xInrun,yInrun)
    axis equal
    grid on
    title('approach: height vs distance')
    xlabel('horizontal distance (m)')
    ylabel('height above starting point (m)')
    hold on
plot(xInrun_raw, yInrun_raw)

dydxInrun= [0 diff(yInrun)./diff(xInrun)];
ddydxInrun=[0 diff(dydxInrun)];
kurvature=ddydxInrun./(1+dydxInrun.^2).^1.5;

% given variables
CdA=0.557*0.8 % m^2 0.557 is seated position from Hoerner book
rho0=1.1839; T0=298.15; Y0=8772.; T=0+T0;

Y=3870; % elevation of Loveland in meters
rho=rho0*exp(-T0*Y/T/Y0)*T0/T
g=9.81
m=60
mu=0.05
eta=(CdA*rho)/2/m

tF=7 ;   % final time in seconds                         
dt=0.001;     % time step
tspan=0:dt:tF;

options = odeset('Events',@takeoff_event2);
% states are [x     v  ]
 state0=     [0     1.2]; %initial condition
%integrate ode's 10
ode=@(t,state) inrun2(t,state);
[ti,statei,tei,yei,iei]=ode45(ode,tspan,state0,options);

tyiei=[tei yei iei]

x=statei(:,1);
v=statei(:,2)*mps2mph;

figure(2)
plot(x,v)
    axis equal    
    grid on
    title('speed vs. distance')
    xlabel('distance (meters)')
    ylabel('speed (miles per hour)')

my function is:

    function [ statedot ] = inrun2(t,state)
    % ODEs for inrun acceleration
    % state vector = [ x  v ]

    global xInrun yInrun dydxInrun kurvature
    global m CdA mu eta g
    x=state(1);% x = horizantal coordinate
    v=state(2);% v = horizantal velocity
    dydx=interp1(xInrun,dydxInrun,x);% finding dydx for any value of x
    k=interp1(xInrun,kurvature,x);% finding curvature for any value of x
    theta=atan(dydx);% angle of slope at any value of x
    statedot=[0 0]';% setting initial value of statedot
    statedot(1)=v*cos(theta);
    N=m*(g*cos(theta)+k*v*v);% normal force
    statedot(2)=-g*sin(theta)-eta*v*v-mu*N*sign(v)/m; 
    txvstatedot=[t x dydx v statedot']
    %pause
    end

    % x = horizantal coordinate
    % v = velocity
    % state[x v]
    % statedot[xdot(horizantal velocity) vdot(acceleration along surface)]
    % N = normal force

[Jan]

Without a proper formatting, the code is not readable. Please follow the "? Help" link to learn how to use the forum's interface.

[Dean]

I'm sorry for making it so messy. I reformatted it, so I hope that makes it easier to see.

[Jan]

@Dean: Thank you very much. You are one of the very few persons (when I remember about 5 since the forum has been started!), who react personally to the suggestion to format the code. Therefore I've voted for your question. See Answers: Code formatting in the forum.

Answer 1

[Jan]

You define:

tF=7 ;   % final time in seconds                         
dt=0.001;     % time step
tspan=0:dt:tF;

This should be sufficient already for a fine resolution of the output. But I do not understand, what "x=53 to x=59" means, and it is nearly impossible to read the unformatted code. What is "x"? When you specify tspan as a [1 x 2] vector and use low relative and absolute tolerances, you get a finer resoltuion of the output, when the stepsize control of the integrator detects large changes in the trajectory. But then the time-steps are not uniformly anymore.

[Dean]

x is the horizontal component of the position- i fixed my question to be a little more readable...sorry about that, it's my first time. Anyway, you'll see from figure 1 that the section that ode45 is skipping over is the transition which, dynamically speaking, is the most interesting part and I want good data from that region. So, I can set tspan as a [1x2] vector but how do I set the solver to use low relative and absolute tolerances? And what is considered a low tolerance?

[Jan]

As usual in Matlab, you find the corresponding instructions, when you read the documentation:

help ode45
doc ode45

Look for "options" set by odeset, parameters RelTol, AbsTol and NormControl. They are used by the step size control: ODE45 calculates two trajectories, one with an order 4 and one with an order 5 Runge-Kutta-scheme. If both trajectories have a relative distance smaller than RelTol and an absolute distance smaller than AbsTol, a step is accepted. If the tolerances are not met, the step size is halfed. The step size for the next step is adjusted according to the distances to keep the total number of steps as small as possible.

You cannot define uniquely, what "low" is. It depends on the magnitudes of the processed data. It matters, if you simulate a flight to mars or a bouncing ball in your room in the units meters. If the limits are chosen too small, the number of calculated steps grows and the accumumated rounding errors will reduce the accuracy of the result. If the tolerances are too high, the large step sizes reduce the accuracy mainly by the large discretization errors.

Therefore choosing the tolerances, deciding if NormControl is helpful and scaling the data (using good units) is hard science (or purely arbitrary from a another point of view). While tiny test problems like predator/prey systems can be solved easily by beginners also using ODE45, it would be a bad idea to simulate a large and complicated systems like a 10'000 degree of freedom complete vehicel simulation naively just by calling ODE45 with the default parameters.

Actually a trajectory determined by an integrator is not a scientific solution already, but it is required to apply small variations of the inputs and parameters to control the sensitivity of the results. E.g. a simulation of a pencil standing on its tip is not meaningful, because the final position is extremely sensitive to the initial conditions. It would be even better not to vary the initial conditions, but to perform this variation in each step and accumulate the sesitivity matrices (also known as "Wronski"-matrix), but this is not implemented in Matlab's ODE45.