How can I plot the z-displacement of my floater correctly ?

This time I have difficulties to let the following code work. It has bugs and I don't know how to solve it. This time I want to plot the displacement of a floater on Faraday waves. Thank you very much for your help !

% Define simulation parameters
rho = 1000; % fluid density [kg/m^3]
mu = 1e-3; % fluid viscosity [Pas]
g = 9.81; % acceleration due to gravity [m/s^2]
L = 1; % length of container [m]
W = 0.5; % width of container [m]
A = 0.1; % amplitude of oscillation [m]
omega = 2*pi; % frequency of oscillation [rad/s]
rho_f = 500; % floater density [kg/m^3]
V = 0.01; % floater volume [m^3]
% Set up grid
Nx = 100;
Ny = 50;
x = linspace(0, L, Nx); % intervalle 0 - L divisé en Nx segments
y = linspace(0, W, Ny);
[X, Y] = meshgrid(x, y);
dx = L/Nx; % grid spacing in x direction
dy = W/Ny; % grid spacing in y direction
% Set up time-stepping
dt = 0.001; % time step [s]
tmax = 10; % maximum time [s]
t = 0:dt:tmax; % time vector
% Set up initial conditions
u = zeros(Ny, Nx); % initial x velocity [m/s]
v = zeros(Ny, Nx); % initial y velocity [m/s]
eta = zeros(Ny, Nx); % initial displacement [m]
% Set up initial position of floaters
x_floater = L/2; % initial x position of floater [m]
y_floater = W/2; % initial y position of floater [m]
% Run simulation
for i = 2:length(t)
    % Compute acceleration at current time step
    [a, b] = acceleration(u, v, eta, rho, mu, g, A, omega, t(i), rho_f, V, x_floater, y_floater);
    % Update velocity and displacement using Euler method
    u = u + adt;
    v = v + bdt;
    eta = eta + udt + vdt;
    % Update position of floater
    x_floater = x_floater + u(round(y_floater/dy), round(x_floater/dx))*dt;
    y_floater = y_floater + v(round(y_floater/dy), round(x_floater/dx))*dt;
    % Check if floater is within bounds of container
    if x_floater < 0 || x_floater > L
        u_floater = -u_floater; % reverse x velocity of floater
    end
    if y_floater < 0 || y_floater > W
        v_floater = -v_floater; % reverse y velocity of floater
    end
end

Unrecognized function or variable 'dx'.

Error in solution>acceleration (line 62)
[eta_x, eta_y] = gradient(eta, dx, dy);

% Visualize results
figure;
surf(X, Y, eta);
xlabel('x');
ylabel('y');
zlabel('displacement');
title('Faraday waves');
hold on;
plot3(x_floater, y_floater, eta(round(y_floater/dy), round(x_floater/dx)), 'r.', 'MarkerSize', 20);
hold off;
% Define function to compute acceleration
function [a, b] = acceleration(u, v, eta, rho, mu, g, A, omega, t, rho_f, V, x_floater, y_floater)
% Compute x and y gradients of displacement
[eta_x, eta_y] = gradient(eta, dx, dy);
% Compute x and y components of velocity at floater position
u_floater = interp2(X, Y, u, x_floater, y_floater);
v_floater = interp2(X, Y, v, x_floater, y_floater);
% Compute x and y components of acceleration at floater position
a_floater = -rho*g*eta_x(round(y_floater/dy), round(x_floater/dx)) - 2*mu*u_floater/V - rho_f*A*omega^2*sin(omega*t);
b_floater = -rho*g*eta_y(round(y_floater/dy), round(x_floater/dx)) - 2*mu*v_floater/V - rho_f*A*omega^2*sin(omega*t);
% Interpolate acceleration at each grid point
a = interp2(X, Y, a_floater, X, Y);
b = interp2(X, Y, b_floater, X, Y);
end

Answer 1

[Walter Roberson]

You need to pass dx and dy into the function acceleration() -- or you need to pass in enough information to be able to calculate the values. You also need to pass in X and Y.

(Provided that you can guarantee that your X and Y passed in to the functions are regular grids, you could calculate dx and dy from the X and Y values.)