Error in ODE arguments

Hello,

I am a newbie at MATLAB and I am trying to get a forcing function to work. I can get it to work without the F(t)/M but once I add it, it says error in ODE arguments:

clear

close all

% System parameters

m1 = 2000; Icg = 2500; % kg

c1 = 3000; c2 = 3000; % kg/s

k1 = 30000; k2 = 30000; % N/m

l1 = 1; l2 = 1.5;

M = [m1 0

0 Icg];

C = [(c1+c2) (c1*l1-c2*l2)

(c1*l1-c2*l2) ((c2*l2^2)+(c1*l1^2))];

K = [(k1+k2) (k1*l1-k2*l2)

(k1*l1-k2*l2) ((k2*l2^2)+(k1*l1^2))];

Br = [k1 k2

k1*l1 -k2*l2];

Brdot = [c1 c2

c1*l1 -c2*l2];

r = [.015 .015]';

r1 = [.086 .086]';

F = @(t) Br*r + Brdot*r1;

% Time grid

t0 = 0; tf = 10; dt = 0.01; t = t0:dt:tf;

% Set initial state and integrate equations of motion

s0 = [1 1 0 0]';

f = @(t,s) [s(3);s(4);(F(t)/M)*-C/M*[s(3) s(4)]'-K/M*[s(1) s(2)]'];

[t,s] = ode45(f,t,s0);

% Plot system motion

figure

subplot(211),plot(t,s(:,1),t,s(:,2),'LineWidth',2)

grid minor

legend('x1','x2')

xlabel('Time (s)')

ylabel('Displacements (m)')

title('System Dynamic Response')

subplot(212),plot(t,s(:,3),t,s(:,4),'LineWidth',2)

grid minor

legend('v1','v2')

xlabel('Time (s)')

ylabel('Velocity (m/s)')

[Jan]

Please remove the pile of blank lines from your code and use the environment for code to improve the readability. In addition we can run your code by copy&paste afterwards.

Can you explain, where the blank lines are coming from? I see them in the forum from time to time but do not understand, how they are created.

Whenever you mention an error in the forum, post a copy of the complete error message. This is a very useful information and you have it already on the screen. It is easier to solve a problem, when the problem is known.

[Alexander Salas]

%   Numerical solution of IVP

%   M*xddot + C*xdot + K*x = F(t)  with x(0) = x0 and xdot(0) = v0

 

% System parameters

m1 = 2000; Icg = 2500; % kg

c1 = 3000; c2 = 3000; % kg/s

k1 = 30000; k2 = 30000; % N/m 
l1 = 1; l2 = 1.5;
M = [m1 0

     0 Icg];

C = [(c1+c2) (c1*l1-c2*l2)

     (c1*l1-c2*l2)    ((c2*l2^2)+(c1*l1^2))];

K = [(k1+k2) (k1*l1-k2*l2)

     (k1*l1-k2*l2)   ((k2*l2^2)+(k1*l1^2))];

Br = [k1 k2
      k1*l1  -k2*l2];
Brdot = [c1 c2
         c1*l1  -c2*l2];
r = [.015 .015]';
r1 = [.086 .086]';

F = @(t) (Br*r) + (Brdot*r1);
 
% Time grid

t0 = 0;  tf = 10; dt = 0.01; t = t0:dt:tf;

% Set initial state and integrate equations of motion

s0 = [1 1 0 0]';
f  = @(t,s) [s(3);s(4);F(t)/M-C/M*[s(3) s(4)]'-K/M*[s(1) s(2)]'];
[t,s]  = ode45(f,t,s0);  

% Plot system motion

figure

subplot(211),plot(t,s(:,1),t,s(:,2),'LineWidth',2)
grid minor
legend('x1','x2')
xlabel('Time (s)')
ylabel('Displacements (m)')
title('System Dynamic Response')
subplot(212),plot(t,s(:,3),t,s(:,4),'LineWidth',2)
grid minor
legend('v1','v2')
xlabel('Time (s)')
ylabel('Velocity (m/s)')

[Alexander Salas]

[Jan]

This is no valid Matlab syntax:

M = [m1 0
    
     0 Icg];

The blank line in the code is not allowed. Therefore I asked you to remove the blank lines.

Answer 1

[Alan Stevens]

I suspect you mean like this (notice the way M divides, using the back-slash):

%   Numerical solution of IVP
%   M*xddot + C*xdot + K*x = F(t)  with x(0) = x0 and xdot(0) = v0
% System parameters

m1 = 2000; Icg = 2500; % kg
c1 = 3000; c2 = 3000; % kg/s
k1 = 30000; k2 = 30000; % N/m 
l1 = 1; l2 = 1.5;
M = [m1 0
     0 Icg];
C = [(c1+c2) (c1*l1-c2*l2)
     (c1*l1-c2*l2)    ((c2*l2^2)+(c1*l1^2))];
K = [(k1+k2) (k1*l1-k2*l2)
     (k1*l1-k2*l2)   ((k2*l2^2)+(k1*l1^2))];
Br = [k1 k2
      k1*l1  -k2*l2];
Brdot = [c1 c2
         c1*l1  -c2*l2];
r = [.015 .015]';
r1 = [.086 .086]';

F = @(t) (Br*r) + (Brdot*r1);
 
% Time grid

t0 = 0;  tf = 10; dt = 0.01; t = t0:dt:tf;

% Set initial state and integrate equations of motion

s0 = [1 1 0 0]';
f  = @(t,s) [s(3);s(4);M\F(t)-M\C*[s(3) s(4)]'-M\K*[s(1) s(2)]']; %%%%%%%%%%%%%%
[t,s]  = ode45(f,t,s0);  

% Plot system motion

figure

subplot(211),plot(t,s(:,1),t,s(:,2),'LineWidth',2)
grid minor
legend('x1','x2')
xlabel('Time (s)')
ylabel('Displacements (m)')
title('System Dynamic Response')
subplot(212),plot(t,s(:,3),t,s(:,4),'LineWidth',2)
grid minor
legend('v1','v2')
xlabel('Time (s)')
ylabel('Velocity (m/s)')

[Jan]

A further simplification:

f  = @(t,s) [s(3); s(4); M \ F(t) - M \ C * s(3:4) - M \ K * s(1:2)];