6 views (last 30 days)

Show older comments

Hi,

I've been reading various posts on getting extra parameters from ODEs and not having much luck implementing them. Concerned that I'm trying something too complicated for my Matlab skills, I've gone back to a basic model to try and understand where I'm going wrong. This has raised a couple of questions.

First, the code, the majority of which is borrowed from others. This is the top half of a simple two degree of freedom mass-spring-damper model (without the plotting and analysis parts):

clear

clc

% define time span

t0 = 0; % s. Input signal start time

t1 = 5.115; % s. Input signal end time. 5.115 gives 1024 points; 61.435 gives 12,288 points

dt = 0.005; % time resolution

tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

% define chirp input

A = 10; % mm. Input signal peak amplitude

f0 = 0.5; % Hz. Input signal start frequency

f1 = 20; % Hz. Input signal end frequency

g = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency

i = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal

idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

% set initial conditions

x0 = [0; 0; 0; 0];

% Sprung mass parameters

ms = 540.5; % kg

ks = 41; % N/mm

cs = 1.5; % Ns/mm

% Unsprung mass parameters

mu = 40; % kg

ku = 350; % N/mm

cu = 0.35; % Ns/mm

% Solve model

[T, x] = ode45(@(t,x) Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu), tvec, x0);

And the ODE function is:

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

i = interp1(tvec, i, t); % this is interpolating i at t

idot = interp1(tvec, idot, t); % this is interpolating idot at t

% as Matlab is using unspecified time steps it needs a value of i for

% each t

% x(1) is sprung mass displacement

% x(2) is unsprung mass displacement

dx(1) = x(3); % sprung mass velocity. This is first column of dx

dx(2) = x(4); % unsprung mass velocity. This is second column of dx, and so on....

dx(3) = 1./ms.*(ks.*(x(2)-x(1)) + cs.*(x(4)-x(3))).*1000; % sprung mass acceleration

dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - (ks.*(x(2)-x(1))) - (cs.*(x(4)-x(3)))).*1000; % unsprung mass acceleration

dx = dx'; % transpose results from horizontal to vertical

end

This runs fine but the first question is why the need to do the transpose at the bottom when I don't see that line in anyone else's examples? Without it, the code throws an error, and I'm aware that the ODE must return the results in columns. However, I'm confused as to why I don't see the transpose in the help pages or here on Answers.

I then adjusted the code to see if I could get an extra parameter out of the ODE - just a simple dummy parameter as an example. I inserted the following in the ODE:

Fs = ks.*(x(2)-x(1)) + cs.*(x(4)-x(3)); % sum of forces on sprung mass

And adjusted two of the ODE statements to accept the new parameter:

dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration

dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - Fs).*1000; % unsprung mass acceleration

I adjusted the function declaration to include the new parameter:

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

I then asked for the parameter from the ODE:

[dx, Fs] = Two_DOF_QCM_Basic_ODE(T, x, i, idot, tvec, ms, ks, cs, mu, ku, cu);

It ran without error but I only got Fs back with a single value in it, rather than a value for each element of T. How do I get Fs back at all times in T?

Thanks, Simon.

Mischa Kim
on 19 Jan 2021

Edited: Mischa Kim
on 19 Jan 2021

Hi Simon, to your first question: When you assign values to dx(1), dx(2), and so on you are creating a row vector. However, you need for the ode45 function to return a column vector, hence you have to transpose. What is frequently done in the ode function is to intialize the vector, e.g.

dx = zeros(4,1); % this is a 4x1 column vector

as a column vector. Assigning values to the dx now preserves the column vector. The other approach (which I use frequently) is to assign dx as:

dx = [x(3); ... % sprung mass velocity. This is first column of dx

x(4); ... % unsprung mass velocity. This is second column of dx, and so on....

1./ms.*(ks.*(x(2)-x(1)) + cs.*(x(4)-x(3))).*1000; ... % sprung mass acceleration

1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - (ks.*(x(2)-x(1))) - (cs.*(x(4)-x(3)))).*1000]; % unsprung mass acceleration

which by design is a column vector. To your second question, since you have an equation for Fs you can simply calculate Fs after solving the ODE, i.e. with the return values from the ode45 call.

Stephen
on 25 Jan 2021

Edited: Stephen
on 25 Jan 2021

Using the OutputFcn is a complex way to get the Fs values.

The simpler approach is to run the ODE solver normally, and then run the function with the values returned by the ODE solver, which automatically calculates Fs values which correspond exactly to the T and x values returned by the ODE solver. This just takes two lines of code, as shown below:

% define time span

t0 = 0; % s. Input signal start time

t1 = 5.115; % s. Input signal end time. 5.115 gives 1024 points; 61.435 gives 12,288 points

dt = 0.005; % time resolution

tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

% define chirp input

A = 10; % mm. Input signal peak amplitude

f0 = 0.5; % Hz. Input signal start frequency

f1 = 20; % Hz. Input signal end frequency

g = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency

i = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal

idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

% set initial conditions

x0 = [0; 0; 0; 0];

% Sprung mass parameters

ms = 540.5; % kg

ks = 41; % N/mm

cs = 1.5; % Ns/mm

% Unsprung mass parameters

mu = 40; % kg

ku = 350; % N/mm

cu = 0.35; % Ns/mm

% Solve model

fun = @(t,x) Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu);

[T, x] = ode45(fun, tvec, x0) % solve

[~,Fs] = cellfun(fun,num2cell(T),num2cell(x,2),'uni',0); % This is all you need...

Fs = cell2mat(Fs) % ... and the Fs values correspond exactly to the T and x values.

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

i = interp1(tvec, i, t); % this is interpolating i at t

idot = interp1(tvec, idot, t); % this is interpolating idot at t

% as Matlab is using unspecified time steps it needs a value of i for

% each t

Fs = ks.*(x(2)-x(1)) + cs.*(x(4)-x(3)); % sum of forces on sprung mass

dx = nan(4,1); % !!! preallocate column output !!!

% x(1) is sprung mass displacement

% x(2) is unsprung mass displacement

dx(1) = x(3); % sprung mass velocity. This is first column of dx

dx(2) = x(4); % unsprung mass velocity. This is second column of dx, and so on....

dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration

dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - Fs).*1000; % unsprung mass acceleration

end

Stephen
on 25 Jan 2021

"...l'm a bit confused about your statement that the outputfcn approach collects all intermediate points..."

Sorry, that was my mistake (since removed from my answer): I confused the common attempt of using persistent variables or similar to store all intermediate calculations, including those backward steps, etc., with the use of outputfcn.

I would have to check the outputFcn documentation, and probably try some examples.

In any case, it is still the more complex approach, if your goal is just to get the Fs values.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!