Can't solve ODE with array instead of constant

Hello there,

I am currently working on a project where I am trying to model a vehicle in my mechanical vibrations class. I have developed a code that works when the external forces are constants, but when I introduce the external forces as arrays, I get some errors. The specific point is in line 10 (the long one). The two arrays are hx and ht (both 1x1001). If I set ht and hx equal to some constant, it solves easily, but once I try with the arrays it won't solve the problem.

Do you have any tips? Thanks in advance

/Kristian

 A = importdata('Realroad.m');
m=5000;
M=75000;
kh=72417751;
ks=1384941;
ch=10000;
cs=100000;
f = @(t,x) [x(3);
    x(4);
    (-ch/m)*x(3)+(ch/m)*x(4)-(cs/m)*x(3)-(kh/m)*x(1)+(kh/m)*x(2)-(ks/m)*x(1)+hx.*(ks/m)+ht.*(cs/m);
    (ch/M)*x(3)-(ch/M)*x(4)+(kh/M)*x(1)-(kh/M)*x(2);];
[t,xa] = ode45(f,[0:0.1:100],[0 0 0 0]);
subplot(3,2,1)
plot(t,xa(:,1))
title('pos')
grid on
xlabel('t')
ylabel('m')
subplot(3,2,2)
plot(t,xa(:,2))
title('pos')
grid on
xlabel('t')
ylabel('m')
subplot(3,2,3)
plot(t,xa(:,3))
title('vel')
grid on
xlabel('t')
ylabel('m/s')
subplot(3,2,4)
plot(t,xa(:,4))
title('vel')
grid on
xlabel('t')
ylabel('m/s')
dy=diff(xa(:,3))./diff(t);
subplot(3,2,5)
plot(t(2:end),dy/9.81)
title('acc')
grid on
xlabel('t')
ylabel('g')
dy=diff(xa(:,4))./diff(t);
subplot(3,2,6)
plot(t(2:end),dy/9.81)
title('acc')
grid on
xlabel('t')
ylabel('g')

end

[John BG]

Hi Kristian

would it be possible for you to attach the data Realroad.m?

It greatly assists readers interested in your question to avoid guessing whether the problem may come from other places like, the missing data.

[Kristian Therkildsen]

Hi John,

Of course - I am currently working on it, so it is not final yet. Ht is therefore not included yet, only hx. The script is a plot of a terrain profile. Thanks in advance :-)

    clear all;close all;
% spatial frequency (n0) cycles per meter
Omega0 = 0.1;
% psd ISO (used for formula 8)
Gd_0 =  32768 * (10^-6);
% waveviness
w = 2; 
% road length
L = 250;
%delta n 
N = 1000;
Omega_L = 0.004;
Omega_U = 1;
delta_n =  1/L; % delta_n = (Omega_U - Omega_L)/(N-1);
% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;
%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);
% calculate amplitude using formula(8) in the article
%Amp = sqrt(2*Gd*delta_n);
%calculate amplitude using simplified formula(9) in the article
k = 6;
Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
%random phases 
Psi = 2*pi*rand(size(Omega)); 
% x abicsa from 0 to L
x = 0:0.25:250;
% road sinal
hx= zeros(size(x));
for i=1:length(x)
    hx(i) = sum( Amp.*cos(2*pi*Omega*x(i) + Psi) );
end
subplot(2,1,1)
plot(x,hx*1000);
xlabel('Distance m');
ylabel('Elevation (mm)');

/Kristian

[Jan]

The code is strange:

A = importdata('Realroad.m');

It is unusual, that an M-file contains data, which can be imported. It is not clear, where ht is coming from. In consequence I cannot run your code. The current description "it won't solve the problem" does not allow to understand, what you observe. Do you get an error message? If so, post a copy of the complete message.

Answer 1

[Jan]

Do not solve an integration for a large set of parameters at once, because this reduces the processing speed and accuracy of the result. Remember, that the step-size controller of the integrator chooses the time steps such, that the effect of local discretization errors and the accumulated rounding error is minimized. Performing the integration with a set of different parameters does not allow to run the integration at the optimal trajectory for each parameter.

To run ODE for a set or different parameters, call it in a loop and vary the parameters. Using an anonymous function for the actual integration is not useful: slower and harder to debug. Use a standard function instead:

function dx = fcn(t, x, ch, m, cs, kh, hx, ht, M)
dx = [x(3); ...
      x(4);
      (-ch/m)*x(3) + (ch/m)*x(4) - (cs/m)*x(3) - (kh/m)*x(1) + ...
      (kh/m)*x(2) - (ks/m)*x(1) + hx.*(ks/m) + ht.*(cs/m); ...
      (ch/M)*x(3) - (ch/M)*x(4) + (kh/M)*x(1) - (kh/M)*x(2)];
end

Now the loop:

for hxi = 1:numel(hx)
  for hti = 1:numel(ht)
    fcnP   = @(t,x) fcn(t, x, ch, m, cs, kh, hx(hxi), ht(hti), M);
    [t,xa] = ode45(fcnP, 0:0.1:100, [0 0 0 0]);
    ...
  end
end

[Kristian Therkildsen]

Hi Jan,

Thanks for your help. I am kind of new to MATLAB, so I still have a lot to learn, but I can feel the progress. I finished the realroad.m file in order to get both hx and ht.

I also edited my previous script because of your advice and tried to assemble my scripts, but even with your tips, I still get an error. Perhaps it's an easy mistake, but I am having trouble solving it. I get the following:

Error: File: newscript.m Line: 75 Column: 1 Function definitions in a script must appear at the end of the file. Move all statements after the "fcn" function definition to before the first local function definition. Line 75 is where hxi starts. I really appreciate the help :-)

/Kristian

   patial frequency (n0) cycles per meter
Omega0 = 0.1;
% psd ISO (used for formula 8)
Gd_0 =  32768 * (10^-6);
% waveviness
w = 2; 
% road length
L = 250;
%delta n 
N = 1000;
Omega_L = 0.004;
Omega_U = 1;
delta_n =  1/L; % delta_n = (Omega_U - Omega_L)/(N-1);
% spatial frequency band
Omega = Omega_L:delta_n:Omega_U;
%PSD of road
Gd = Gd_0.*(Omega./Omega0).^(-w);
% calculate amplitude using formula(8) in the article
%Amp = sqrt(2*Gd*delta_n);
%calculate amplitude using simplified formula(9) in the article
k = 6;
Amp = sqrt(delta_n) * (2^k) * (10^-3) * (Omega0./Omega);
%random phases 
Psi = 2*pi*rand(size(Omega)); 
% x abicsa from 0 to L
x = 0:0.25:250;
% road sinal
hx= zeros(size(x));
for i=1:length(x)
    hx(i) = sum( Amp.*cos(2*pi*Omega*x(i) + Psi) );
end
subplot(2,1,1)
plot(x,hx*1000);
xlabel('Distance m');
ylabel('Elevation (mm)');
max(Amp)
%Velocity 10 km/h = 2.778m/s
v=2.778;
t=x./v;
ht=diff(hx)./diff(t);
subplot(2,1,2)
plot(t(2:end),ht*1000)
xlabel('time t');
ylabel('velocity mm/s');
%Second script
m=5000;
M=75000;
kh=72417751;
ks=1384941;
ch=10000;
cs=100000;
function dx = fcn(t, x, ch, m, cs, kh, hx, ht, M)
dx = [x(3); ...
      x(4);
      (-ch/m)*x(3) + (ch/m)*x(4) - (cs/m)*x(3) - (kh/m)*x(1) + (kh/m)*x(2) - (ks/m)*x(1) + hx.*(ks/m) + ht.*(cs/m);
      (ch/M)*x(3) - (ch/M)*x(4) + (kh/M)*x(1) - (kh/M)*x(2)];
end
for hxi = 1:numel(hx)
    for hti = 1:numel(ht)
    [t,xa] = ode45(fcnP, 0:0.1:100, [0 0 0 0]);
    fcnP   = @(t,x) fcn(t, x, ch, m, cs, kh, hx(hxi), ht(hti), M);
    end
end
subplot(4,2,1)
plot(t,xa(:,1))
title('pos')
grid on
xlabel('t')
ylabel('m','Rotation',0)
subplot(4,2,2)
plot(t,xa(:,2))
title('pos')
grid on
xlabel('t')
ylabel('m','Rotation',0)
subplot(4,2,3)
plot(t,xa(:,3))
title('vel')
grid on
xlabel('t')
ylabel('m/s','Rotation',0)
subplot(4,2,4)
plot(t,xa(:,4))
title('vel')
grid on
xlabel('t')
ylabel('m/s','Rotation',0)
dy=diff(xa(:,3))./diff(t);
subplot(4,2,5)
plot(t(2:end),dy/9.81)
title('acc')
grid on
xlabel('t')
ylabel('g','Rotation',0)
dy=diff(xa(:,4))./diff(t);
subplot(4,2,6)
plot(t(2:end),dy/9.81)
title('acc')
grid on
xlabel('t')
ylabel('g','Rotation',0)