Quarter Car Model with Haversine Bump Input, Errors using ode45

Hey all,

I am trying to simulate a quarter car model, using a user defined Euler function (Also having trouble with ode45). But no matter what I try I keep getting errors:

Here is my main code:

%Speed Bump Conditons %This values while be variable in the GUI, the fixed
%numbers are for testing
Height = 0.1;               %Height of the bump         [m] 
Length = 0.6;                %Length of the bump         [m]
Velocity = 40;              % Velocity of the car     [m/s]

%Initial conditions, assumed to be zero 
vs = 0 ;                        %Sprung mass velocity      [m/s]
vt = 0;                         %Unsprung mass velocity    [m/s]
ys = 0 ;                        %Initial sprung mass height [m]
yt = 0;                         %Intial unsprung mass height [m]

%Values used for Euler's Method
tspan = [0 3];                      %Span of time for simulation [seconds]
x = [vs, vt, ys, yt];               %Initial conditions 
h = 0.00001;                          % Defining Step Size
t = tspan(1):h:tspan(2);


%%%Haversine Bump Input
Bump =  Height*(sin(pi*Velocity*t/Length)).^2;

for k = 1:length(t)
    if t(k) > T
        Bump(k) = 0; %Setting bump height values less than 0 for
        % time greater than T, the time it takes to travel over bump
    end %ending if statement
end %Ending for loop

% %Defining the anonymous function
% [dy] = ODECalc(Variables,Bump');
k1 = Variables(1);
m1 = Variables(2);
k2 = Variables(3);
c2 = Variables(4);
m2 = Variables(5);

dy = @(t,y) [(1/m2) * ((-k2*y(2)) + (k2*y(4)) - (c2*y(1)) + (c2*y(3)));
                y(1);
              (1/m1) * (((-k2-k1)*y(4)) - (c2*y(3)) + (c2*y(1)) + (k2*y(2)) + (k1*Bump'));
                y(3)];
%[t2,x2] = ode45(dy,t,x)    %Doesn't work
%Solving the ODE using Eulers method
[t2,x2] = odeEuler(dy,tspan,h,x);

The code for the Eulers function is:

function [x,y] = odeEuler(f,range,h,y0)
%This function solves the ODE using Euler's method. The function allows a
%system of first order ODEs to be solved. 

%Function Inputs:
%   f : handle to a MATLAB function, 
%   range : range of independited variables of which the simulation will be
%       performed.
%   h : step size for the range
%   y0 : inital values for the ODE
%Function Outputs:
%   x : row vector of independent variables of the approximate solution
%   from Euler's Method
%   y : The soltuion array that corresponds with the x vector

%Begin function code.

% Initialize the output solution arrays.
x = (range(1):h:range(2)); %Intializing x vector, from begining of range, to end of range with step size h
y = y0'; %Initializing y vector, as initial condtion, y0, must be transposed inorder to perform calculations

% Compute x an y using Euler's method. 
for i = 1:length(x)-1 %start for loop
    y(:,i+1) = y(:,i) + f(x(:,i),y(:,i))*h; %Algorithm to calulate the next y values
end %Ending for loop



y = y';     % output y as a series of column vectors.
end %End of function 

The error I get using this configuration is:

Arrays have incompatible sizes for this operation.

Error in odeEuler (line 24)

y(:,i+1) = y(:,i) + f(x(:,i),y(:,i))*h; %Algorithm to calulate the next y values

Error in Calculations (line 58)

[t2,x2] = odeEuler(dy,tspan,h,x);

If I define the equation bump at an anonymous function, I get results however it is not correct as the input continues, whereas I want it to stop after the define length of the bump.

Any help would be greatly appreciated.

Answer 1

[Walter Roberson]

Bump is a row vector. k1 is a scalar. k1*Bump' is a column vector. So your dy is producing a column vector that is longer than your tspan, and it is doing that for every call, for every individual time. The resulting vector is the wrong size for the rest of the computation

[Walter Roberson]

I suggest that you make Bump into a function that calculates the Havershine from time passed in and multiplies by (time<T) to zero out for larger times

[Walter Roberson]

But remember that your function must have continuous first and second derivatives, so suddenly zeroing out Bump might be a problem. A way to work with that would be to run the ode to time T and then run it again for all remaining time.

[Brendan]

Thank you for the reply Walter,

Can you elaborate a little on how I would do that? I think I get the jist of what you are saying to do but do not know how I would implement it. I orginally tried having it as a function, but still got errors.

My original bump function essentially just how I did the calculation above. As expected I recieved the same errors.

function [Bump, T] = BumpCalc(Height, Length, Velocity,t)
%This function calculates the bump input dependant on the inputted length
%height and velocity the user inputs. 
%The time it takes to travel over the bump is calculated 
%The value for the bump is outputted, as well as the time it takes to
%travel over the bump


T = Length/Velocity;        %Time to travel over the bump [s]


%%%%%%%%%%%%%%%%%%%%% Havershine Bump Equation %%%%%%%%%%%%%%%%%%%%%%%%%%
Bump = Height*(sin(pi*Velocity*t/Length)).^2;

%%%%%%%%%%% Setting variables equal to zero, past time T %%%%%%%%%%%%%%%%%
for k = 1:length(t)
    if t(k) > T
        Bump(k) = 0; %Setting bump height values less than 0 for
        % time greater than T, the time it takes to travel over bump
    end %ending if statement
end %Ending for loop
end

[Walter Roberson]

T = Length/Velocity; 
Bump = @(t) Height*sinpi(t./T).^2 .* (t <= T);

dy = @(t,y) [(1/m2) * ((-k2*y(2)) + (k2*y(4)) - (c2*y(1)) + (c2*y(3)));
                y(1);
              (1/m1) * (((-k2-k1)*y(4)) - (c2*y(3)) + (c2*y(1)) + (k2*y(2)) + (k1*Bump(t)));
                y(3)];

Question: are you sure you only want a half cycle, sin(pi*fraction), not a full cycle sin(2*pi*fraction) ?

[Brendan]

Thank you so much Walter, I was pulling my hair out trying to figue this out. To answer you question I do indeed want it to be a half cycle, the input simulates the car driving over 1 bump. Once again thank you.