CODEWITHSUNDEEP
matlabodenumerical-integration
freddie_k
freddie_k

Reputation: 13

MATLAB - passing a sinusoidal forcing function to ode45

I'm new to Matlab and am really struggling even to get to grips with the basics.

I've got a function, myspring, that solves position and velocity of a mass/spring system with damping and a driving force. I can specify values for the spring stiffness (k), damping coefficient (c), and mass (m), in the command window prior to running ode45. What I am unable to do is to define a forcing function (even something simple like g = sin(t)) and use that as the forcing function, rather than having it written into the myspring function.

Can anyone help? Here's my function:

function pdot = myspring(t,p,c,k,m)

w = sqrt(k/m);

g = sin(t);  % This is the forcing function

pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = g - c*p(2) - (w^2)*p(1);

end

and how I'm using it in the command window:

>> k = 2; c = 2; m = 4;
>> tspan = linspace(0,10,100);
>> x0 = [1 0];
>> [t,x] = ode45(@(t,p)myspring(t,p,c,k,m),tspan,x0);

That works, but what I want should look something like this (I imagine):

function pdot = myspring(t,p,c,k,m,g)

w = sqrt(k/m);

pdot = zeros(size(p));
pdot(1) = p(2);
pdot(2) = g - c*p(2) - (w^2)*p(1);

end

Then

g = sin(t);
[t,x] = ode45(@(t,p)myspring(t,p,c,k,m,g),tspan,x0);

But what I get is this

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in myspring (line 7)
pdot(2) = g - c*p(2) - (w^2)*p(1);

Error in @(t,p)myspring(t,p,c,k,m,g)

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Horchler, thank you for the reply. I can do as you suggested and it works. I am now faced with another problem that I hope you could advise me on.

I have a script that calculates the force on a structure due to wave interaction using Morison's equation. I gave it an arbitrary time span to begin with. I would like to use the force F that script calculates as the input driving force to myspring. Here is my Morison script:

H = 3.69;           % Wave height (m)
A = H/2;            % Wave amplitude (m)
Tw = 9.87;           % Wave period (s)
omega = (2*pi)/Tw;  % Angular frequency (rad/s)
lambda = 128.02;    % Wavelength (m)
k = (2*pi)/lambda;  % Wavenumber (1/m)
dw = 25;             % Water depth (m)
Cm = 2;             % Mass coefficient (-)
Cd = 0.7;           % Drag coefficient (-)
rho = 1025;         % Density of water (kg/m^3)
D = 3;              % Diameter of structure (m)

x = 0;              % Fix the position at x = 0 for simplicity


t = linspace(0,6*pi,n);   % Define the vector t with n time steps
eta = A*cos(k*x - omega*t); % Create the surface elevation vector with size equal to t

F_M = zeros(1,n); % Initiate an inertia force vector with same dimension as t
F_D = zeros(1,n); % Initiate a drag force vector with same dimension as t
F = zeros(1,n); % Initiate a total force vector with same dimension as t

fun_inertia = @(z)cosh(k*(z+dw));  % Define the inertia function to be integrated
fun_drag = @(z)(cosh(k*(z+dw)).*abs(cosh(k*(z+dw))));   % Define the drag function to be integrated

for i = 1:n
    F_D(i) = abs(((H*pi)/Tw) * (1/sinh(k*dw)) * cos(k*x - omega*t(i))) * ((H*pi)/Tw) * (1/sinh(k*dw)) * cos(k*x - omega*t(i)) * integral(fun_drag,-dw,eta(i));
    F_M(i) = (Cm*rho*pi*(D^2)/4) * ((2*H*pi^2)/(Tw^2)) * (1/(sin(k*dw))) * sin(k*x - omega*t(i)) * integral(fun_inertia,-dw,eta(i));
    F(i) = F_D(i) + F_M(i);
end

Any further advice would be much appreciated.

Upvotes: 1

Views: 1714

Answers (1)

horchler
horchler

Reputation: 18484

You can't pre-calculate your forcing function. It depends on time, which ode45 determines. You need to define g as a function and pass a handle to it into your integration function:

...
g = @(t)sin(t);
[t,x] = ode45(@(t,p)myspring(t,p,c,k,m,g),tspan,x0);

And then call it I n your integration function, passing in the current time:

...
pdot(2) = g(t) - c*p(2) - (w^2)*p(1);
...

Upvotes: 2

Related Questions