Solve equations of motion containing matrices and partial derivatives

Hello,

I am having an issue running a simulation for my equation of motion.

I am calculating a potential/electric field over some values of x and y (a matrix). I want to use it afterwards to compute my equations of motion which are the following:

with

I have tried many things but still no result. For example, I tried to run it using while and for loops but no success. I am also not very good at using Runge Kutta 4 so I could not make it work for this complex case.

Would anybody be able to help?

Here is my code for the electric field (first part, note that a_p is d_p in the code):

% Initialize variables
nbE = 4;
lambda = 4;
g = 1;
e = 35e-3;
w = 1;
f = 50;
omega = 2 * pi * f;
V0 = 1000;
% grid_size = 100;
x_len = 8;
y_len = 2.5;
grid_size = 200;
dx = x_len / grid_size;
dy = y_len / grid_size;
x_min = 0; x_max = 8;
y_min = -0.0; y_max = 2.5;
yLoc = 0;
d = 0;
x = linspace(x_min, x_max, grid_size);
y = linspace(y_min, y_max, grid_size);
K1 = 10;


% Discretization parameters
dx = 0.08;
dy = 0.08;
x = 0:dx:8;
y = e:dy:2.5;

% Initialize matrices
[V, Ex, Ey] = deal(zeros(length(y), length(x), length(t)));
A = zeros(length(y));

t = 0;
% Iterate over grid
for it = 1:length(t)
    for ix = 1:length(x)
        for iy = 1:length(y)
            if(y(iy) < 0)
                A(iy) = 1/2 * V0 * K1 * exp((2*pi*y(iy))/lambda);
            elseif(y(iy) == 0)
                A(iy) = 1/2 * V0 * K1 * exp((-2*pi*(y(iy)+e))/lambda);
            else
                A(iy) = 1/2 * V0 * K1 * exp((-2*pi*y(iy))/lambda);
            end
            V(iy, ix, it) = A(iy) * cos((2*pi*x(ix)/lambda) - omega*t(it)) + A(iy) * cos((2*pi*x(ix)/lambda) + omega*t(it));
            Ex(iy, ix, it) = (2*pi/lambda) * A(iy) * sin((2*pi*x(ix)/lambda) - omega*t(it)) + (2*pi/lambda) * A(iy) * sin((2*pi*x(ix)/lambda) + omega*t(it));
            Ey(iy, ix, it) = (2*pi/lambda) * A(iy) * cos((2*pi*x(ix)/lambda) - omega*t(it)) + (2*pi/lambda) * A(iy) * cos((2*pi*x(ix)/lambda) + omega*t(it));
        end
    end
end

and for my equations:

epsilon_0 = 8.854*10^(-12); % Free space permittivity
gM = 1.62*10^3; % g_z
d_p = 3*10^(-2);
rho_bulk = 1.50*10^(-3);
sigma_s = 2.64*10^(-11);
V_p = 4/3 * pi * (d_p/2)^3;
m_p = V_p * rho_bulk;
beta = 0; % Lagging phase in radians
lambda = 4 * 10^(-3); 
V0 = 1000; % Electric potential amplitude in Volts
Ec = 3*10^3;
d = 1.25*10^(-3);
mu = 1.85*10^(-5);
z0 = 3*10^(-7);
r = d_p;
ws = 3.8*10^(-8);

epsilon_glass = 4;
mu_g = 1.85 * 10^(-5);

% Compute H
hama = 12 * pi * z0^2 * ws;

%% Simulation Parameters
% Compute partial derivatives of Ex and Ey
[dEx_dy, dEx_dx] = gradient(Ex);
[dEy_dy, dEy_dx] = gradient(Ey);

t_start = 0;
t_end = 250; % In ms
dt = 1;  % Time step (in ms)
dt(1) = dt;
y0 = 4.3e-3;  % Initial x position
f = 5;
tspan = t_start:dt:t_end; % Time span (in ms)
% Lengthspan in y axis
ye = linspace(0, 20, 4);


qs = 4*pi*(d_p/2)^2 * sigma_s;
qp = 0.1 * qs;

% Compute constants
T = 0;
C = qp/m_p;
K = ((4*pi*(d_p/2)^3)/m_p) * epsilon_0 * ((epsilon_rp - 1)/(epsilon_rp + 2));
M = qp^2/(16*pi*epsilon_0*m_p);

% Initial conditions
y0 = 4.3;
y = [];
z = [];
vy = [];
vz = [];

y(1) = 4.3;
z(1) = 35e-3;
vy(1) = 0;
vz(1) = 0;
ay(1) = 0;
az(1) = 0;

% Discretization parameters
dy = 0.08; % Spatial step in x-direction (mm)
dz = 0.08; % Spatial step in y-direction (mm)
y_range = 0:dy:8; % y-dimension from 0 to 8 mm
z_range = 35e-3:dz:4; % (zLoc-d) z-dimension from 0 to 4 mm

% Preallocate solution
sol = zeros(4,numel(tspan));
it = 1;

% Initial conditions
initialConditions = [y(1); vy(1); z(1); vz(1)];

% Time points
tspan = [0 t_end];

% Split solution
y_sol = y(:,1);
vy_sol = y(:,2);
z_sol = y(:,3); 
vz_sol = y(:,4);

a = 1;

while t < t_end
    t
    % Find the index with the closest value in the electric field ranges
    [~, idx_y] = min(abs(y_range - y(it)));
    [~, idx_z] = min(abs(z_range - z(it)));

    % Compute y double dot
    ay(it+1) = K*(Ex(idx_z,idx_y)*dEx(idx_z,idx_y) + Ey(idx_z,idx_y)*dEy(idx_z,idx_y)) ...
             +  C*Ex(idx_z,idx_y);
    % Compute z double dot
    az(it+1) = K*(Ex(idx_z,idx_y)*dEx(idx_z,idx_y) + Ey(idx_z,idx_y)*dEy(idx_z,idx_y)) ...
             + C*Ey(idx_z,idx_y) - hama/(6*m_p) * ((r * (d_p/2))/(z(1)^2*(r + (d_p/2))) + (d_p/2)/(z(1) + r)^2) ...
             - M/z(it) - gM;

    % Integrate using RK4
    % Initial conditions
    y0 = y(it);
    vy0 = vy(it);
    z0 = z(it);
    vz0 = vz(it);

    % RK4 slopes
    k1y = vy0;
    k1vy = ay(it);

    k1z = vz0;
    k1vz = az(it);

    k2y = dt(it) * (vy(it) + k1vy*1/2);
    k2vy = dt(it) * ay(it) + k1z*(1/2);

    k2z = vz0 + k1vz*dt/2;
    k2vz = dt(it) * (az(it) + (1/2)*k1z);
    k3y = dt(it) *(k2vy*1/2);
    k3vy = dt(it) * (ay(it) + k2z*1/2);

    k3z = dt(it) *(vz(it) + (1/2)*k2vz);
    k3vz = dt(it) * (az(it) + k2z*1/2);

    k4y = dt(it) * (vy(it) + k3vy);
    k4vy = dt(it) * (ay(it) + k3y);

    k4z = dt(it) * (vz(it) + k3vz);
    k4vz = dt(it) * (az(it) + k3z);

    % Update next state
    y(it+1) = y(it,1) + (dt(it)/6)*(k1y + 2*k2y(1) + 2*k3y(1) + k4y(1));
    vy(it+1) = vy(it,1) + (dt(it)/6)*(k1vy(1) + 2*k2vy(1) + 2*k3vy(1) + k4vy(1));

    z(it+1) = z(it,1) + (dt(it)/6)*(k1z(1) + 2*k2z(1) + 2*k3z(1) + k4z(1));  
    vz(it+1) = vz(it,1) + (dt(it)/6)*(k1vz(1) + 2*k2vz(1) + 2*k3vz(1) + k4vz(1));

    % [y, vy, z, vz] = ODEIntegration(t_start, t_end, dt(i), [y(1) z(1)], [vy(1) vz(1)], [ay az], 1e-8);

    if (y(it+1) - y(it))/y(it+1) < 0.01 && ((z(it+1) - z(it))/z(it+1) < 0.01)
        dt(it+1) = dt(1);
    else
        dt(it+1) = dt(it)/10;
    end

    t = t + dt(it)
end

Thank you!

Answer 1

[Sulaymon Eshkabilov]

Your posted code has a few critical issues to be addressed before executing it. (1) epsilon_rp is undefined; (2) y(1) =4.3; but the cose is requesting y to be of a size 1 by 4; (3) dEx is neither defined nor computed; (4) According your shown formulations, dEy and dEz need to be computed; etc.