Search code examples
matlabphysics

How can I add a function of time requested previously as input in ODE 45? [MATLAB]

I'm having this problem: Let's suposse that a user wants to add an arbitrary Force (a function of time) to a spring mass damper system. In order to do that, we request the force as an input, i.e. F = cos(t). Well, I noticed that if we add a function of time in the code where the F is placed (inside ODE), then the solutions for the movement equations are correct. But i can't find a way to achieve this goal. I mean, to request the force as input so that ODE changes everytime i put a new force. I've tried almost everything like converting the input F to a function, etc.

I hope you could help me. Here is the problem in the code:

F   = input('Insert force F(t)..............[N] = ','s'); 

u0  = [x0,v0];                                            % INITIAL CONDITIONS AS A VECTOR

ODE = @(t,u) [u(2) ; F - (c/m)*u(2) - (k/m)*u(1)];        % FORCED SPRING MASS DAMPER SYSTEM

[t,u] = ode45(ODE, [0 tf], u0);

The fix code results:

m   = input('Insert mass of the system..........................[kg] = ');
c   = input('Insert viscosity coefficient........................[ ] = ');
k   = input('Insert elasticity constant..........................[m] = ');
x0  = input('Insert initial position.............................[m] = ');
v0  = input('Insert initial velocity...........................[m/s] = ');
tf  = input('Insert time of complete movement....................[s] = ');
F   = input('Insert external force as F(t).......................[N] = ','s');

F = str2func( ['@(t)' F] );                    % CONVERT STRING TO FUNCTION

%=========================================================== ODE 45 =======
u0 = [x0,v0];                              % INITIAL CONDITIONS AS A VECTOR

ODE  = @(t,u) [u(2) ; F(t) - (c/m)*u(2) - (k/m)*u(1)];
[t,u] = ode45(ODE, [0 tf], u0);

u(:,3) = F(t) - (c/m)*u(:,2) - (k/m)*u(:,1);                 % ACCELERATION

%============================================================= PLOT =======
v = {'x','dx/dt','d2x/dt2'};                         % CHAR VECTOR FOR PLOT

for i = 1:3                                                % PLOT WITH LOOP
    subplot(3,1,i)
    plot(t,u(:,i));
    xlabel('t');
    ylabel(v(i));
    grid on;
end

Code just for personal uses. Not public because input allows malicious commands.

Thanks to @Wolfie and @Ander Biguri.

Solution

  • I'm going to take the liberty of simplifying your problem down to its main parts.

    1. You have some char array from your user input which defines a function in terms of t. I'm going to entirely ignore all the reasons that blindly evaluating user-input functions is dangerous, but you shouldn't.

    2. You want to combine this function with an existing function in terms of t to produce a single function

    We can do this like so:

    F = 'cos(t)'; % char version of function, e.g. from the input command
y = @(t) t;   % some base function in terms of t

F = str2func( ['@(t)' F] ); % Convert function char to anonymous function

yF = @(t) y(t) + F(t); % combine the functions

    You can test this with a quick plot

    tt = 0:0.01:10; % axis data
figure; hold on;
plot( tt, y(tt),  'displayname', 'y(t)'  );
plot( tt, yF(tt), 'displayname', 'yF(t)' );
legend('show');

    plot

    In your actual example you have ODE in place of my example f, so assuming F is converted to an anonymous function as shown above, you would get

    ODE = @(t,u) [u(2) ; F(t) - (c/m)*u(2) - (k/m)*u(1)];