Solve equations of motion containing matrices and partial derivatives

7 次查看（过去 30 天）
Styrbjorn
Styrbjorn 2024-1-30
回答： Sulaymon Eshkabilov 2024-1-31
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!

请先登录，再回答此问题。

回答（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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by