MATLAB

Question

我是 Matlab 的新手，我什至很难掌握基础知识。

我有一个函数 myspring，它可以求解具有阻尼和驱动力的 mass/spring 系统的位置和速度。我可以在 window 之前的命令中指定 spring 刚度 (k)、阻尼系数 (c) 和质量 (m) 的值运行ode45。我无法做的是定义一个强制函数（甚至像 g = sin(t) 这样简单的东西）并使用 that 作为强制函数，而不是将其写入 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

以及我如何在命令中使用它 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);

可行，但我想要的应该是这样的（我想）：

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

但我得到的是这个

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

霍希勒，谢谢你的回复。我可以按照你的建议去做，而且它有效。我现在面临另一个问题，希望你能给我一些建议。

我有一个脚本可以使用莫里森方程计算由于波浪相互作用而作用在结构上的力。一开始我给了它一个任意的时间跨度。我想使用脚本计算的力 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);

然后在您的集成函数中调用它，传入当前时间：

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

MATLAB - 将正弦强制函数传递给 ode45

MATLAB - passing a sinusoidal forcing function to ode45

ode

numerical-integration