Plotting the Evolution of Spatially Localized Initial Conditions for Discrete Klein-Gordon Equation

4 次查看（过去 30 天）
Athanasios Paraskevopoulos
回答： Athanasios Paraskevopoulos 2024-6-10
Ran in:
Ιn continuation of the study with which I have been fascinated Numerical Simulation of a Damped, Driven Nonlinear Wave System with Spatially Extended Initial Conditions
Now we examine the role of damping in the structure of the branches of equilibrium points. In this study, we consider localized initial conditions of the form:
where is the distance between the nodes, and is the amplitude. We assume zero initial velocity
We also note that if we consider an infinite chain, the initial conditions $(1)$ satisfy the Dirichlet boundary conditions only asymptotically. However, since the smallest half-length of the chain has been set to $ L/2 = 100$, the error is of the order of $ 10^{-44} $, which does not affect numerical computations.
The figure shows the dynamics of the initial condition (1) for parameters: and . In the first image of Figure, we observe the profile of the initial condition.
The subsequent images show the dynamics of the system for different values of the damping force . For , the convergence takes place at the equilibrium point which is the basic state.
I want to create the red and the blue plots
Now I have tried but as you can see the result are different. What should I change?
% Parameters
L = 200; % Length of the system
K = 99; % Number of spatial points
omega_d = 1; % Characteristic frequency
beta = 1; % Nonlinearity parameter
delta_values = [0.1, 0.05, 0.01]; % Damping coefficients
% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+1); % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
% Time parameters
dt = 1; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector
% Values of amplitude 'a' to iterate over
a_values = [2]; % Initial amplitude for the initial condition plot
% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
U = Y(1:N);
Udot = Y(N+1:end);
Uddot = zeros(size(U));
% Laplacian (discrete second derivative)
for k = 2:N-1
Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) ;
end
% System of equations
dUdt = Udot;
dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);
% Pack derivatives
dYdt = [dUdt; dUdotdt];
end
% Create a figure for subplots
figure;
% Loop through each value of 'delta' and generate the plot
for i = 1:length(delta_values)
delta = delta_values(i);
deltaScaled = h * delta;
a = a_values(1);
% Initial conditions
U0 = a * sech(-L/2 + n*h); % Initial displacement
U0(1) = 0; % Boundary condition at n = 0
U0(end) = 0; % Boundary condition at n = K+1
Udot0 = zeros(size(U0)); % Initial velocity
% Pack initial conditions
Y0 = [U0, Udot0];
% Solve ODE
opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
[t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);
% Extract solutions
U = Y(:, 1:N);
Udot = Y(:, N+1:end);
% Plot final displacement profile
subplot(2, 2, i+1);
plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title(['$t=3000$, $\delta=', num2str(delta), '$'], 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-2 2]);
grid on;
end
% Initial plot
subplot(2, 2, 1);
a_init = 2; % Example initial amplitude for the initial condition plot
U0_init = a_init * sech(-L/2 + n*h); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;
% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed

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

采纳的回答

Athanasios Paraskevopoulos
Ran in:
% Parameters
L = 200; % Length of the system
K = 99; % Number of spatial points
omega_d = 1; % Characteristic frequency
beta = 1; % Nonlinearity parameter
delta_values = [0.1, 0.05, 0.01]; % Damping coefficients
% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+2); % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
% Time parameters
dt = 1; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector
% Values of amplitude 'a' to iterate over
a_values = [1]; % Initial amplitude for the initial condition plot
% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
U = Y(1:N);
Udot = Y(N+1:end);
Uddot = zeros(size(U));
% Laplacian (discrete second derivative)
for k = 2:N-1
Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) ;
end
% System of equations
dUdt = Udot;
dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);
% Pack derivatives
dYdt = [dUdt; dUdotdt];
end
% Create a figure for subplots
figure;
% Loop through each value of 'delta' and generate the plot
for i = 1:length(delta_values)
delta = delta_values(i);
deltaScaled = h * delta;
a = a_values(1);
% Initial conditions
U0 = a * sech( n*h); % Initial displacement
U0(1) = 0; % Boundary condition at n = 0
U0(end) = 0; % Boundary condition at n = K+1
Udot0 = zeros(size(U0)); % Initial velocity
% Pack initial conditions
Y0 = [U0, Udot0];
% Solve ODE
opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
[t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);
% Extract solutions
U = Y(:, 1:N);
Udot = Y(:, N+1:end);
% Plot final displacement profile
subplot(2, 2, i+1);
plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title(['$t=3000$, $\delta=', num2str(delta), '$'], 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-2 2]);
grid on;
end
% Initial plot
subplot(2, 2, 1);
a_init = 1; % Example initial amplitude for the initial condition plot
U0_init = a_init * sech( n*h); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;
% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Line Plots 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by