# Solve equations of motion containing matrices and partial derivatives

7 views (last 30 days)
Styrbjorn on 30 Jan 2024
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
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!

Sulaymon Eshkabilov on 31 Jan 2024
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.

### Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2023a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!