MATLAB - передача синусоидальной форсирующей функции ode45

Question

Я новичок в Matlab, и я действительно борюсь даже за то, чтобы справиться с основами.

У меня есть функция, myspring, которая решает положение и скорость системы массы/пружины с демпфированием и движущей силой. Я могу указать значения жесткости пружины (k), коэффициента демпфирования (c) и массы (m) в командном окне до запуска ode45. То, что я не могу сделать, это определить функцию принуждения (даже что-то простое, например g = sin(t)), и использовать , чтобы в качестве функции принуждения вместо того, чтобы записать ее в функцию myspring.

Может ли кто-нибудь помочь? Вот моя функция:

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 

и как я использую его в окне командной строки:

>> 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); 

Это работает, но то, что я хочу, должен выглядеть примерно так (я представляю себе):

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 

Тогда

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

Bu t, что я получаю, это

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, спасибо за ответ. Я могу сделать то, что вы предложили, и это работает. Теперь я столкнулся с другой проблемой, на которую, надеюсь, вы могли бы мне посоветовать.

У меня есть сценарий, который вычисляет силу на структуре из-за взаимодействия волн с использованием уравнения Морисона. Я дал ему произвольный промежуток времени для начала. Я хотел бы использовать силу F, которую скрипт вычисляет как движущую силу ввода на myspring. Вот мой Morison сценарий:

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 

Любые дальнейшие рекомендации были бы оценены.

Answer 1

Вы не можете предварительно вычислить свою функцию принуждения. Это зависит от времени, которое определяет ode45. Вам нужно определить g как функцию и передать a handle to it в вашу функцию интеграции:

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

А потом называть его I п вашу интеграционную функцию, передавая в текущее время:

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