From 3a4edcb31206e0f8abaa23f3285d4dbbf43b027b Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 1 Apr 2025 12:13:54 +0200 Subject: [PATCH] New test and example scripts. --- Estimations/CustomCylindricalCutOff.m | 77 +++++++++ Estimations/EstimateCostFunctionMLoop.m | 22 ++- Estimations/MinimizeEnergyFunctionalViaCGD.m | 165 ++++++++++++++++++- Estimations/SplitStepFourierMethod.m | 91 ++++++++++ 4 files changed, 346 insertions(+), 9 deletions(-) create mode 100644 Estimations/CustomCylindricalCutOff.m create mode 100644 Estimations/SplitStepFourierMethod.m diff --git a/Estimations/CustomCylindricalCutOff.m b/Estimations/CustomCylindricalCutOff.m new file mode 100644 index 0000000..589a6de --- /dev/null +++ b/Estimations/CustomCylindricalCutOff.m @@ -0,0 +1,77 @@ +% Define the size of the cube using NumberOfGridPoints +NumberOfGridPoints = [32, 32, 16]; % [X, Y, Z] dimensions of the cube +Dimensions = [30, 30, 15]; % Real-world dimensions for the grid + +% Generate a random 3D cube based on the NumberOfGridPoints +cube = rand(NumberOfGridPoints(1), NumberOfGridPoints(2), NumberOfGridPoints(3)); + +% Define the real-world distance range for each axis (X, Y, Z) +x_range = linspace(-Dimensions(1)/2, Dimensions(1)/2, NumberOfGridPoints(1)); % 32 points along X mapped to [-15, 15] +y_range = linspace(-Dimensions(2)/2, Dimensions(2)/2, NumberOfGridPoints(2)); % 32 points along Y mapped to [-15, 15] +z_range = linspace(-Dimensions(3)/2, Dimensions(3)/2, NumberOfGridPoints(3)); % 16 points along Z mapped to [-7.5, 7.5] + +% Create a 3D grid of coordinates for the cube using meshgrid +[x, y, z] = meshgrid(x_range, y_range, z_range); + +% Define the cylinder's radius and height (in real-world units) +radius = 10; % Radius in real-world units (e.g., meters) +height = 4; % Height in real-world units (e.g., meters) + +% Get the center of the cube (in real-world coordinates) +center_x = mean(x_range); % Center of x_range +center_y = mean(y_range); % Center of y_range +center_z = mean(z_range); % Center of z_range + +% Define the cylindrical mask +cylinder_mask = ((x - center_x).^2 + (y - center_y).^2 <= radius^2) & ... + (abs(z - center_z) <= height / 2); + +% Set everything inside the cylinder to 1 and outside to 0 +mask = double(cylinder_mask); + +% Multiply the mask by the random cube to apply the cut-off +masked_cube = cube .* mask; + +% Visualize the cube and the cylindrical mask using isosurfaces +figure(1); +clf +figure_width = 1200; +set(gcf, 'Position', [100, 100, figure_width, round(figure_width / 1.618)]); +set(gca, 'FontSize', 16, 'Box', 'On', 'Linewidth', 2); + +% Set the axis limits to be the same for both subplots +x_limits = [min(x_range), max(x_range)]; +y_limits = [min(y_range), max(y_range)]; +z_limits = [min(z_range), max(z_range)]; + +% Visualize the random cube +subplot(1, 2, 1); +isosurface(x, y, z, cube, 0.5); % Isosurface for the random cube +axis([x_limits, y_limits, z_limits]); +title('Original 3D Numerical Grid', 'FontSize', 16); +xlabel('X (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +ylabel('Y (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +zlabel('Z (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +grid on; +axis equal; +view(3); +camlight; +lighting gouraud; + +% Visualize the masked cube +subplot(1, 2, 2); +isosurface(x, y, z, masked_cube, 0.5); % Isosurface for the masked cube +axis([x_limits, y_limits, z_limits]); +title('Grid with Custom Cylindrical Cut-off', 'FontSize', 16); +% Add subtitle with radius and height with adjusted position +text('Units', 'normalized', 'Position', [0.5, -0.2], 'String', ... + ['Radius = ', num2str(radius), ' \mum, Height = ', num2str(height), ' \mum'], ... + 'Interpreter', 'tex', 'FontSize', 14, 'FontWeight', 'normal', 'HorizontalAlignment', 'center'); +xlabel('X (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +ylabel('Y (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +zlabel('Z (\mum)', 'Interpreter', 'tex', 'FontSize', 16); +grid on; +axis equal; +view(3); +camlight; +lighting gouraud; \ No newline at end of file diff --git a/Estimations/EstimateCostFunctionMLoop.m b/Estimations/EstimateCostFunctionMLoop.m index 9d8747f..f0ddb44 100644 --- a/Estimations/EstimateCostFunctionMLoop.m +++ b/Estimations/EstimateCostFunctionMLoop.m @@ -1,6 +1,6 @@ %% % Parameters -alpha = 0.5; % Example value for alpha +alpha = 8; % Example value for alpha N = linspace(1e4, 1e6, 100); % Atom number range OD = linspace(0, 3, 100); % Optical depth range N1 = 100; % Example parameter for N1 controlling transition steepness @@ -37,13 +37,10 @@ title(['|C(X)| = OD^3 * N^{\alpha - 9/5} * Smoothed Heaviside Function with Alph % title('Cost Function: OD^3 * N^{\alpha - 9/5} with Smoothed Heaviside Function'); %% % Parameters -N = linspace(1e3, 1e6, 100); % Atom number range -OD = linspace(0, 3, 100); % Optical depth range +N = linspace(1e3, 1e6, 100); % Atom number range +OD = linspace(0, 3, 100); % Optical depth range N1 = 1e5; % Example parameter for N1 controlling transition steepness -% Smoothed Heaviside function as given in the appendix - - % Create figure and initial plot fig = figure(2); clf @@ -51,17 +48,26 @@ set(gcf,'Position',[50 50 950 750]) set(gca,'FontSize',16,'Box','On','Linewidth',2); alpha = 0.5; % Initial value for alpha +function cf = getSmoothedHeaviside(x, N1) + if x > 0 + cf = (2 ./ (1 + exp(N1 ./ x))); + else + cf = 0; + end +end + % Nested function to calculate and update the cost function plot function update_plot(alpha_value, N, N1, OD) - SmoothedHeaviside = @(x) (x > 0) .* (2 ./ (1 + exp(N1 ./ x))); + % Calculate the cost function for the given alpha CostFunction = zeros(length(N), length(OD)); for i = 1:length(N) for j = 1:length(OD) - H_N = SmoothedHeaviside(N(i)); % Smoothed Heaviside function + H_N = getSmoothedHeaviside(N(i), N1); % Smoothed Heaviside function CostFunction(i, j) = H_N * OD(j)^3 * N(i)^(alpha_value - 9/5); end end + % Create a meshgrid for plotting [OD_grid, N_grid] = meshgrid(OD, N); % Plot the cost function as a surface plot diff --git a/Estimations/MinimizeEnergyFunctionalViaCGD.m b/Estimations/MinimizeEnergyFunctionalViaCGD.m index 51d42d6..113b1b8 100644 --- a/Estimations/MinimizeEnergyFunctionalViaCGD.m +++ b/Estimations/MinimizeEnergyFunctionalViaCGD.m @@ -1,3 +1,165 @@ +%% Direction solution to the unconstrained optimization problem with the Conjugate Gradient technique +clc; +clear; +format long; + +% Define the function (now accepts a 3D array) +f = @(X) sum(X(:).^4) - 4*sum(X(:).^3) + 6*sum(X(:).^2) + 5; + +% Initial Guess: Now we use a 3D array +x_o = rand(3, 3, 3); % Random 3D array as the initial guess +fprintf('Initial Objective Function Value: %.6f\n\n', f(x_o)); + +% Convergence Criteria: +Epsilon = 1E-5; + +% Gradient Calculation (numerical) +J = compute_gradient(f, x_o); % Gradient at initial point +S = -J; % Initial search direction + +% Iteration Counter: +i = 1; + +% Minimization +while norm(S(:)) > Epsilon % Halting Criterion + % Step Size Optimization (Line Search) + lambda = optimize_step_size(f, x_o, S); + + % Update the solution + x_o_new = x_o + lambda * S; + + % Update the gradient and search direction + J_old = J; + J = compute_gradient(f, x_o_new); + S = update_search_direction(S, J, J_old); + + % Update for next iteration + x_o = x_o_new; + i = i + 1; +end + +% Output +fprintf('Number of Iterations for Convergence: %d\n\n', i); +fprintf('Objective Function Minimum Value: %.6f\n\n', f(x_o)); + + +%% Visualize the optimization + +clc; +clear; +format long; + +% Define the function (now accepts a 3D array) +f = @(X) sum(X(:).^4) - 4*sum(X(:).^3) + 6*sum(X(:).^2) + 5; + +% Initial Guess: Now we use a 3D array +x_o = rand(3, 3, 3); % Random 3D array as the initial guess +fprintf('Initial Objective Function Value: %.6f\n\n', f(x_o)); + +% Convergence Criteria: +Epsilon = 1E-5; + +% Gradient Calculation (numerical) +J = compute_gradient(f, x_o); % Gradient at initial point +S = -J; % Initial search direction + +% Iteration Counter: +i = 1; + +% Prepare for visualization +fig = figure; +clf +set(gcf,'Position',[50 50 950 750]) +set(gca,'FontSize',16,'Box','On','Linewidth',2); +hold on; + +% Create a mesh grid for visualization +[X1, X2] = meshgrid(-2:0.1:2, -2:0.1:2); % Adjust domain based on your function +Z = arrayfun(@(x1, x2) f([x1, x2, 0]), X1, X2); % Example 2D slice: f([X1, X2, 0]) + +% Plot the surface +surf(X1, X2, Z); +xlabel('X1'); +ylabel('X2'); +zlabel('f(X)'); +title('Conjugate Gradient Optimization Process'); +colorbar; +scatter3(x_o(1), x_o(2), f(x_o), 'r', 'filled'); % Initial point + +% Set the 3D view +view(75, 30); % Example 3D view: azimuth 45°, elevation 30° + +% Minimization with visualization +while norm(S(:)) > Epsilon % Halting Criterion + % Step Size Optimization (Line Search) + lambda = optimize_step_size(f, x_o, S); + + % Update the solution + x_o_new = x_o + lambda * S; + + % Update the gradient and search direction + J_old = J; + J = compute_gradient(f, x_o_new); + S = update_search_direction(S, J, J_old); + + % Update for next iteration + x_o = x_o_new; + i = i + 1; + + % Visualization Update + scatter3(x_o(1), x_o(2), f(x_o), 'g', 'filled'); % Current point + drawnow; % Update the plot +end + +% Output +fprintf('Number of Iterations for Convergence: %d\n\n', i); +fprintf('Objective Function Minimum Value: %.6f\n\n', f(x_o)); + +%% +% Numerical Gradient Calculation using the finite differences method for 3D +function J = compute_gradient(f, X) + epsilon = 1E-5; % Step size for numerical differentiation + J = zeros(size(X)); % Gradient with the same size as the input + + % Loop over all elements of the 3D array + for i = 1:numel(X) + X1 = X; + X2 = X; + % Perturb the element X(i) in both directions + X1(i) = X1(i) + epsilon; + X2(i) = X2(i) - epsilon; + % Central difference formula for gradient + J(i) = (f(X1) - f(X2)) / (2 * epsilon); + end +end + +% Line Search (Step Size Optimization) +function lambda = optimize_step_size(f, X, S) + alpha = 1; % Initial step size + rho = 0.5; % Step size reduction factor + c = 1E-4; % Armijo condition constant + max_iter = 100; % Max iterations for backtracking + + for k = 1:max_iter + % Check if the Armijo condition is satisfied + grad = compute_gradient(f, X); % Compute gradient + if f(X + alpha * S) <= f(X) + c * alpha * (S(:)' * grad(:)) + break; + else + alpha = rho * alpha; % Reduce the step size + end + end + lambda = alpha; % Return the optimized step size +end + +% Update Search Direction (Polak-Ribiere method) +function S_new = update_search_direction(S, J_new, J_old) + beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); + S_new = -J_new + beta * S; +end + + +%{ function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V,Observ) format long; @@ -168,4 +330,5 @@ function S_new = update_search_direction(S, J_new, J_old) % (Polak-Ribiere method) beta = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2)); S_new = -J_new + beta * S; -end \ No newline at end of file +end +%} \ No newline at end of file diff --git a/Estimations/SplitStepFourierMethod.m b/Estimations/SplitStepFourierMethod.m new file mode 100644 index 0000000..04eeeed --- /dev/null +++ b/Estimations/SplitStepFourierMethod.m @@ -0,0 +1,91 @@ +% Constants +hbar = 1.0545718e-34; % Reduced Planck constant (Joule second) +m_e = 9.10938356e-31; % Mass of electron (kg) +wavelength = 5e-9; % de Broglie wavelength (meters) +k0 = 2 * pi / wavelength; % Wave number + +% Grid parameters +Lx = 1e-6; % Simulation domain size (meters) +Nx = 1024; % Number of grid points +dx = Lx / Nx; % Grid spacing +x = linspace(-Lx / 2, Lx / 2, Nx); % Spatial grid + +% Time parameters +dt = 3e-15; % Time step (seconds) +Nt = 1001; % Number of time steps + +% Wave packet parameters +sigma_x = 60e-9; % Width of the Gaussian packet (meters) +x0 = -2e-7; % Initial position of the packet (meters) +kx0 = k0; % Initial wave vector + +% Barrier parameters +barrier_width = 4e-8; % Width of the square barrier (meters) +barrier_height = 0.8 * hbar^2 * k0^2 / (2 * m_e); % Height of the barrier (Joules) + +% Initial Gaussian wave packet +psi0 = exp(-(x - x0).^2 / (2 * sigma_x^2)) .* exp(1i * kx0 * x); +psi0 = psi0 / sqrt(sum(abs(psi0).^2) * dx); % Normalization + +% Potential energy (Square barrier in the center) +V = zeros(size(x)); +V(abs(x) < barrier_width / 2) = barrier_height; + +% Fourier space components +kx = (2 * pi / Lx) * (-Nx/2:Nx/2-1); % Frequency domain representation +kx = fftshift(kx); % Shift zero frequency component to center +K2 = kx.^2; % Squared wave numbers + +% Split-step Fourier method +psi = psi0; % Initialize wave function +transmission_prob = 0; % Initialize transmission probability + +figure(1); % Create a figure for visualization +clf +set(gcf,'Position',[100 100 950 750]) +set(gca,'FontSize',16,'Box','On','Linewidth',2); +for t = 1:Nt + % (a) 1/2 Evolution for the potential energy in real space + psi = psi .* exp(-1i * V * dt / (2 * hbar)); + + % (b) Forward transform to Fourier space + psi_k = fft(psi); + + % (c) Full evolution for the kinetic energy in Fourier space + psi_k = psi_k .* exp(-1i * hbar * K2 * dt / (2 * m_e)); + + % (d) Inverse Fourier transform back to real space + psi = ifft(psi_k); + + % (e) 1/2 Evolution for the potential energy in real space + psi = psi .* exp(-1i * V * dt / (2 * hbar)); + + % Calculate transmission probability + Rho = abs(psi).^2; % Probability density + transmission_prob = sum(Rho(x > barrier_width / 2)) * dx; + + % Visualization of wave packet evolution + if mod(t, 50) == 0 + clf; % Clear figure before drawing + plot(x * 1e9, Rho / max(Rho), 'r', 'LineWidth', 2); % Plot normalized density + hold on; + plot(x * 1e9, V / max(V), 'k', 'LineWidth', 2); % Plot barrier potential + title(sprintf('Time step: %d\nTransmission Probability: %.3f', t, transmission_prob),'FontSize',16); + xlabel('x (nm)', 'FontSize', 16); + ylabel('Probability Density', 'FontSize', 16); + grid on + drawnow; % Update the figure dynamically + end +end + +% Final wave packet distribution +figure(2); +clf +set(gcf,'Position',[100 100 950 750]) +set(gca,'FontSize',16,'Box','On','Linewidth',2); +plot(x * 1e9, abs(psi).^2, 'b', 'LineWidth', 2); +xlabel('x (nm)','FontSize',16); +ylabel('Probability Density','FontSize',16); +title('Final Wave Packet Distribution','FontSize',16); +grid on +%% \ No newline at end of file