Vectorize - C->Matlab / Heat equation

[Koubaros]

Vectorize - C-->Matlab / Heat equation

I want to 'translate' some programs I had in C for MATLAB but by sigtly optimizing the code for MATLAB use. I am new in matlab. So, I had the following code:

(1) for K = 1:dT:M

(2) for I = 1:N+1
dU(I) = (U(I+1)-2*U(I)+U(I-1))/dX^2;
end

(3) for I = 1:N+1
U(I) = U(I)+dU(I)*dT
end

end

This actually is a simple implementation of an explicit scheme for the numerical solution of the heat equation. Now I managed to do this and eliminate the (2) and (3) loop:

(1) for K = 1:dT:M

(2) dU(1:N+1) = (U([2:end 1])-2*U+U([end 1:end-1]))/dX^2;
U(I+1)-->U([2:end 1])
U(I-1)-->U([end 1:end-1])

(3) U(1:N+1) = U(1:N+1)+dU(1:N+1)*dT

end

I could really use your help in order to eliminate the (1) time loop, in a way that the values of U(1:N+1) are stored for each time step dT or, in case M is very large this would take a lot of memory and time, get 'snapshots' every a couple of timesteps which I could define explicitly.

 

[Aaron Bex]

This can be done in MATLAB using a vectorized approach. The following code eliminates the (1) and (2) loops.

U(1:N+1) = U(1:N+1); % Initialize U
for K = 1:dT:M % Loop through each time step
dU(1:N+1) = (U([2:end 1])-2*U+U([end 1:end-1]))/dX^2; % Calculate dU
U(1:N+1) = U(1:N+1)+dU(1:N+1)*dT; % Update U

% Store 'snapshots' of U at every a couple of timesteps
if mod(K,snapshot_timestep) == 0
U_snapshots(:, K/snapshot_timestep+1) = U(1:N+1);
end
end

 

[tacman]

In order to vectorize this code in MATLAB, we can use the built-in functions such as diff and circshift to eliminate the for loops and make the code more efficient. The code can be rewritten as follows:

% Define parameters
dT = 0.01; % time step
dX = 0.1; % spatial step
N = 100; % number of spatial points
M = 1000; % number of time steps

% Initialize U and dU vectors
U = zeros(N+1,1); % initial condition
dU = zeros(N+1,1);

% Time loop
for K = 1:dT:M
% Calculate dU using diff function
dU = diff(U([end 1:end-1])) - 2*U + diff(U([2:end 1]));
dU = dU/dX^2;

% Update U using circshift function
U = U + circshift(dU,1)*dT; % circshift shifts the vector by 1 index

% Save snapshot every 10 time steps
if mod(K,10) == 0
U_snapshots(:,K/10) = U;
end
end

% Plot snapshots
plot(U_snapshots);

In this code, the diff function is used to calculate the difference between neighboring elements in the U vector, eliminating the need for the for loop. Similarly, the circshift function is used to shift the dU vector by 1 index, eliminating the need for the second for loop. Finally, the if statement is used to save snapshots of the solution at every 10 time steps, reducing the memory and time needed to store all the values of U.