Here’s a basic for your MATLAB code:(m-file)

% Parameters

L = 1; % Length of the domain

N = 200; % Number of spatial points (increase for better resolution)

dx = L / (N-1); % Spatial step size

x = linspace(0, L, N); % Spatial grid

nu = 1e-4; % Viscosity

dt = 0.0001; % Time step size (reduce for stability)

T = 1; % Total time

Re = 10000; % Reynolds number

% Initial condition

u = sin(pi*x);

% Time integration using RK4

for t = 0:dt:T

k1 = dt * burgers_rhs(u, nu, dx);

k2 = dt * burgers_rhs(u + 0.5*k1, nu, dx);

k3 = dt * burgers_rhs(u + 0.5*k2, nu, dx);

k4 = dt * burgers_rhs(u + k3, nu, dx);

u = u + (k1 + 2*k2 + 2*k3 + k4) / 6;

end

% Plot the result

figure;

plot(x, u);

xlabel('x');

ylabel('u');

title('Solution of Burgers'' Equation');

% Compute the Fourier transform

U_hat = fft(u);

% Compute the energy spectrum

E = abs(U_hat).^2 / N;

% Plot the energy spectrum

k = (0:N-1) * (2*pi/L);

figure;

loglog(k, E);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

% Function to compute the RHS of Burgers' equation

function rhs = burgers_rhs(u, nu, dx)

N = length(u);

rhs = zeros(size(u));

for i = 2:N-1

rhs(i) = -u(i) * (u(i+1) - u(i-1)) / (2*dx) + nu * (u(i+1) - 2*u(i) + u(i-1)) / dx^2;

end

% Boundary conditions

rhs(1) = 0;

rhs(N) = 0;

end

Thank you very much. But first, I need to non-dimensionalize the equation. I have defined the code as follows, but the results don't seem logical. Can you tell me the reason?

clc

clear

close all

%% Parameter Problem

N_space = 16384; % 32768

x = linspace(0, 1, N_space);

Deltax = x(2) - x(1);

% CFL = 0.5; % convergence condition by Courant–Friedrichs–Lewy

% Deltat = Deltax * CFL;

Deltat = 0.04 * 1e-6;

Tmax = 2;

M_time = Tmax / Deltat; % Cte

t = linspace(0, Tmax, M_time);

Re = 10000;

% k (i)

%% Frist create a zero matrix for the velocity so that the new velocity is calculated and then replaced

% U = zeros(N_space, M_time);

U = zeros(N_space, 1); % Update at 2024-07-13

%% Initial Condition

U(:, 1) = sin(2*pi*x); %%% sin pi*x

%% Boundary Condition

% U(1, :) = 0;

% U(end, :) = 0;

U(1, :) = 0; % Update at 2024-07-13

U(end, :) = 0; % Update at 2024-07-13

%% Initialize figure for plotting

figure(1);

hold on;

tic

%% Implementation rk4

% Spatial Discretization

% Finite Differnce

% dU/dt = 1/Re * Ui+1 - 2*Ui + Ui-1/deltax^2 - 0.5 * Ui+1^2 - Ui-1^2/2*deltax

% f = 1/Re * Ui+1 - 2*Ui + Ui-1/deltax^2 - 0.5 * Ui+1^2 - Ui-1^2/2*deltax

k1 = zeros(N_space, 1);

k2 = zeros(N_space, 1);

k3 = zeros(N_space, 1);

k4 = zeros(N_space, 1);

for n = 1 : M_time

% Un = U(:, n);

Un = U;

for i = 2 : N_space - 1

k1(i) = Deltat * ((-(Un(i+1)^2 - Un(i-1)^2) / (4*Deltax)) + ((1/Re) * (Un(i+1) - 2*Un(i) + Un(i-1)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k2(i) = Deltat * (-(((Un(i+1) + k1(i+1)/2).^2 - (Un(i-1) + k1(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k1(i+1)/2) - 2*(Un(i) + k1(i)/2) + (Un(i-1) + k1(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k3(i) = Deltat * (-(((Un(i+1) + k2(i+1)/2).^2 - (Un(i-1) + k2(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k2(i+1)/2) - 2*(Un(i) + k2(i)/2) + (Un(i-1) + k2(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k4(i) = Deltat * (-(((Un(i+1) + k3(i+1)).^2 - (Un(i-1) + k3(i-1)).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k3(i+1)) - 2*(Un(i) + k3(i)) + (Un(i-1) + k3(i-1))) / (Deltax^2)));

end

for i = 2 : N_space - 1

U(i) = Un(i) + ((1/6)*(k1(i) + 2*k2(i) + 2*k3(i) + k4(i))); % Update at 2024-07-13

end

%% Plot of Result : Plotting U vs x at different times

if n == 1 || mod(n, round(M_time/4)) == 0

plot(x, U, 'DisplayName', ['t = ' num2str(t(n), '%.2f')]);

hold on;

end

%

end

%% Draw Plot

% Drawing velocity plot according to location in 5 different time steps

%% Add labels and legend

xlabel('Space (x)');

ylabel('U (x , t)');

title('Burgers Equation Solution at different times');

legend show;

grid on;

hold off;

%% Calculate energy spectrum

U_hat = fft(U);

E = 1/2*(abs(U_hat)).^2;

k = (0:(N_space/2-1));

% U_hat = fft(U,N_space);

% E = 1/2*((real(U_hat).^2 .* imag(U_hat).^2);

% k = (0:(N_space/2-1));

%% Plot energy spectrum

figure(2);

loglog(k, E(1:(N_space/2)), 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

%% Calculate energy spectrum (Normalization)

N = length(U);

Uhat = fft(U);

Uhat = Uhat / N;

Energy = 1/2*(abs(Uhat)).^2;

Energy = Energy(1:ceil(N/2));

dx = x(2) - x(1);

K = (0:(N/2-1)) / (N * dx);

Energy = Energy / (N * dx);

%% Plot energy spectrum

figure(3);

loglog(K, Energy, 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

SimulationTime = toc

It is possible that there are several reasons why the results are not logical:

Ensure that the boundary conditions are appropriate for the problem you are trying to solve. In the current code, the boundary conditions are zero, which might not be suitable for all cases.

The time step you chose is very small (0.04 * 1e-6). It might be better to increase the time step slightly while ensuring it meets the CFL condition for stability.

The number of spatial points is very large (16384). You can try reducing the number of spatial points to see if that improves the results.

Ensure that the initial conditions are appropriate. In the current code, the initial conditions are a sine function. You can try different initial conditions to see if that improves the results.

Ensure that the function `burgers_rhs` calculates the derivatives correctly. You can add some corrections to improve accuracy.

Try this code:

clc

clear

close all

%% Parameter Problem

N_space = 1024; % Reduced number of spatial points

x = linspace(0, 1, N_space);

Deltax = x(2) - x(1);

Deltat = 0.001; % Increased time step size

Tmax = 2;

M_time = Tmax / Deltat;

t = linspace(0, Tmax, M_time);

Re = 10000;

%% Initialize velocity matrix

U = zeros(N_space, 1);

%% Initial Condition

U(:, 1) = sin(2*pi*x);

%% Boundary Condition

U(1, :) = 0;

U(end, :) = 0;

%% Initialize figure for plotting

figure(1);

hold on;

tic

%% Implementation rk4

k1 = zeros(N_space, 1);

k2 = zeros(N_space, 1);

k3 = zeros(N_space, 1);

k4 = zeros(N_space, 1);

for n = 1 : M_time

Un = U;

for i = 2 : N_space - 1

k1(i) = Deltat * ((-(Un(i+1)^2 - Un(i-1)^2) / (4*Deltax)) + ((1/Re) * (Un(i+1) - 2*Un(i) + Un(i-1)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k2(i) = Deltat * (-(((Un(i+1) + k1(i+1)/2).^2 - (Un(i-1) + k1(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k1(i+1)/2) - 2*(Un(i) + k1(i)/2) + (Un(i-1) + k1(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k3(i) = Deltat * (-(((Un(i+1) + k2(i+1)/2).^2 - (Un(i-1) + k2(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k2(i+1)/2) - 2*(Un(i) + k2(i)/2) + (Un(i-1) + k2(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k4(i) = Deltat * (-(((Un(i+1) + k3(i+1)).^2 - (Un(i-1) + k3(i-1)).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k3(i+1)) - 2*(Un(i) + k3(i)) + (Un(i-1) + k3(i-1))) / (Deltax^2)));

end

for i = 2 : N_space - 1

U(i) = Un(i) + ((1/6)*(k1(i) + 2*k2(i) + 2*k3(i) + k4(i)));

end

%% Plot of Result : Plotting U vs x at different times

if n == 1 || mod(n, round(M_time/4)) == 0

plot(x, U, 'DisplayName', ['t = ' num2str(t(n), '%.2f')]);

hold on;

end

end

%% Draw Plot

xlabel('Space (x)');

ylabel('U (x , t)');

title('Burgers Equation Solution at different times');

legend show;

grid on;

hold off;

%% Calculate energy spectrum

U_hat = fft(U);

E = 1/2*(abs(U_hat)).^2;

k = (0:(N_space/2-1));

%% Plot energy spectrum

figure(2);

loglog(k, E(1:(N_space/2)), 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

%% Calculate energy spectrum (Normalization)

N = length(U);

Uhat = fft(U);

Uhat = Uhat / N;

Energy = 1/2*(abs(Uhat)).^2;

Energy = Energy(1:ceil(N/2));

dx = x(2) - x(1);

K = (0:(N/2-1)) / (N * dx);

Energy = Energy / (N * dx);

%% Plot energy spectrum

figure(3);

loglog(K, Energy, 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

SimulationTime = toc

Abdulkarim Almuhammad

I've sent the results of this code. The energy spectrum plot is correct, but in the velocity plot, the sinusoidal lines should not go above the plot; they should be exactly on the plot, not outside of it.

Another point is that the velocity peak at a certain second should be around 0.8 and should not dissipate too quickly.

Can you help me?

Try this code:

clc

clear

close all

%% Parameter Problem

N_space = 2048; % Increased number of spatial points for better resolution

x = linspace(0, 1, N_space);

Deltax = x(2) - x(1);

Deltat = 0.00001; % Reduced time step size for better stability

Tmax = 2;

M_time = Tmax / Deltat;

t = linspace(0, Tmax, M_time);

Re = 10000;

%% Initialize velocity matrix

U = zeros(N_space, 1);

%% Initial Condition

U(:, 1) = 0.8 * sin(2*pi*x); % Scale initial condition to have peak around 0.8

%% Boundary Condition

U(1, :) = 0;

U(end, :) = 0;

%% Initialize figure for plotting

figure(1);

hold on;

tic

%% Implementation rk4

k1 = zeros(N_space, 1);

k2 = zeros(N_space, 1);

k3 = zeros(N_space, 1);

k4 = zeros(N_space, 1);

for n = 1 : M_time

Un = U;

for i = 2 : N_space - 1

k1(i) = Deltat * ((-(Un(i+1)^2 - Un(i-1)^2) / (4*Deltax)) + ((1/Re) * (Un(i+1) - 2*Un(i) + Un(i-1)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k2(i) = Deltat * (-(((Un(i+1) + k1(i+1)/2).^2 - (Un(i-1) + k1(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k1(i+1)/2) - 2*(Un(i) + k1(i)/2) + (Un(i-1) + k1(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k3(i) = Deltat * (-(((Un(i+1) + k2(i+1)/2).^2 - (Un(i-1) + k2(i-1)/2).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k2(i+1)/2) - 2*(Un(i) + k2(i)/2) + (Un(i-1) + k2(i-1)/2)) / (Deltax^2)));

end

for i = 2 : N_space - 1

k4(i) = Deltat * (-(((Un(i+1) + k3(i+1)).^2 - (Un(i-1) + k3(i-1)).^2) / (4*Deltax)) + ((1/Re) * ...

((Un(i+1) + k3(i+1)) - 2*(Un(i) + k3(i)) + (Un(i-1) + k3(i-1))) / (Deltax^2)));

end

for i = 2 : N_space - 1

U(i) = Un(i) + ((1/6)*(k1(i) + 2*k2(i) + 2*k3(i) + k4(i)));

end

%% Plot of Result : Plotting U vs x at different times

if n == 1 || mod(n, round(M_time/4)) == 0

plot(x, U, 'DisplayName', ['t = ' num2str(t(n), '%.2f')]);

hold on;

end

end

%% Draw Plot

xlabel('Space (x)');

ylabel('U (x , t)');

title('Burgers Equation Solution at different times');

legend show;

grid on;

hold off;

%% Calculate energy spectrum

U_hat = fft(U);

E = 1/2*(abs(U_hat)).^2;

k = (0:(N_space/2-1));

%% Plot energy spectrum

figure(2);

loglog(k, E(1:(N_space/2)), 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

%% Calculate energy spectrum (Normalization)

N = length(U);

Uhat = fft(U);

Uhat = Uhat / N;

Energy = 1/2*(abs(Uhat)).^2;

Energy = Energy(1:ceil(N/2));

dx = x(2) - x(1);

K = (0:(N/2-1)) / (N * dx);

Energy = Energy / (N * dx);

%% Plot energy spectrum

figure(3);

loglog(K, Energy, 'LineWidth', 1);

xlabel('Wavenumber k');

ylabel('Energy E(k)');

title('Energy Spectrum');

grid on;

hold off;

SimulationTime = toc

No it's not correct. Abdulkarim Almuhammad