First commit, adding files

Author: coryhoi

FiniteDifferenceMethods/.gitkeep

@@ -0,0 +1,51 @@
+% Define the parameters
+L = 1;             % Length of the rod
+T = 1;             % Total time
+N = 100;           % Number of spatial grid points
+alpha = 0.02;      % Thermal diffusivity constant
+
+for k=1:5
+    % Refine spacial resolution
+    if k > 1
+        N = N*2;
+    end
+
+    dx = L/N;
+    dt = 0.5*dx*dx/alpha;
+
+    M = int32(T/dt);
+
+    x = linspace(0,L,N+1);
+
+    % Initialize the temperature matrix
+    u = zeros(N+1,M+1);
+
+    % Set the initial condition
+    u(:,1) = sin(pi*x);
+    
+    time = 0.0;
+    toatlError(:) = 0;
+
+    for j = 1:M
+        for i = 2:N
+            u(i,j+1) = u(i,j) + alpha*dt/dx^2 *(u(i+1,j) - 2*u(i,j) + u(i-1,j));
+        end
+    
+        time = time + dt;
+        exact = sin(pi*x)*exp(-alpha*(pi)^2*time);
+    
+        Error = exact' - u(:,j);
+        totalError(j) = sum(Error);
+        
+    end
+    % Error
+    E(k) = norm(totalError,inf);
+    DX(k) = dx;
+end
+
+%%
+figure
+plot(DX,E,'o-')
+xlabel("Delta x")
+ylabel("Error")
+

FiniteDifferenceMethods/convergenceExample.m

@@ -0,0 +1,48 @@
+% Define the parameters
+L = 2*pi;          % Length
+T = 3;             % Total time
+N = 100;           % Number of spatial grid points
+a = 1;             % Wave speed
+
+% Spatial discretization
+dx = L/N;
+
+% Stability condition
+dt = abs(dx/a)*0.5;
+
+% Number of iterations
+M = int32(T/dt);
+
+% Spatial grid
+x = linspace(0,L,N+1);
+
+% Initialize the matrix
+u = zeros(N+1,1);
+
+% Set the initial condition
+u(:,1) = sin(x);
+%%
+
+u_old = u;
+u_older = u;
+
+% Apply the finite difference method
+for i = 1:M
+    % Write your finite difference scheme here
+    % Don't forget about the boundary conditions
+    % And carrying three time levels of the solution
+
+
+
+
+
+    
+end
+
+figure
+plot(x,u)
+xlabel("Position")
+ylabel("U")
+ylim([-1 1])
+xlim([0 2*pi])
+set(gcf,'Position',[400 400 800 380])

FiniteDifferenceMethods/ellipticEquations.mlx

@@ -0,0 +1,94 @@
+% Points in the domain
+nx = 50;
+ny = 50;
+
+% Initializing solution
+U_init = zeros(ny,nx);
+
+% Boundary conditions
+U_init(:,1) = 0;
+U_init(:,nx) = 0;
+U_init(1,:) = 0;
+U_init(ny,:) = 0;
+
+% Corners
+U_init(ny,nx) = 0.5*(U_init(ny-1,nx) + U_init(ny,nx-1));
+U_init(ny,1) = 0.5*(U_init(ny-1,1) + U_init(ny,2));
+U_init(1,nx) = 0.5*(U_init(2,nx) + U_init(1,nx-1));
+U_init(1,1) = 0.5*(U_init(1,2) + U_init(2,1));
+
+U_original = U_init;
+U_new = U_init;
+
+% Optimization parameter
+omega = 2/(1 + (pi/nx));
+
+% Initialize error and iteration count
+maxError = 1;
+iterations = 0;
+
+% Running 3 cases,1 for each method
+numCases = 0;
+
+% Start big loop
+while numCases < 3
+    while maxError > 10^-3
+    
+        for j=2:ny-1
+            for i=2:nx-1
+                % Jacobi method
+                if numCases == 0
+                    U_new(j,i) = 0.25*(U_init(j+1,i) + U_init(j-1,i) + U_init(j,i+1) + U_init(j,i-1));
+                end
+
+                % Gauss-Seidel
+                if numCases == 1
+                    U_new(j,i) = 0.25*(U_new(j+1,i) + U_new(j-1,i) + U_new(j,i+1) + U_new(j,i-1));
+                end
+
+                % SOR
+                if numCases == 2
+                    U_new(j,i) = (1-omega)*U_init(j,i) + (0.25*omega)*(U_new(j+1,i) + U_new(j-1,i) + U_new(j,i+1) + U_new(j,i-1));
+                end
+            end
+        end
+        % Calculate the error
+        
+        % Stopping criteria 1
+        error = U_init - U_new;
+        maxError = max(max(abs(error)));
+
+        % Stopping criteria 2
+        % error = 0;
+        % for j=2:ny-1
+        %     for i=2:nx-1
+        %         error = error + 0.25*(U_new(j+1,i) + U_new(j-1,i) + U_new(j,i+1) + U_new(j,i-1)) - U_new(j,i);
+        %     end
+        % end
+        % maxError = abs(error);
+        
+        % Set U for next iteration
+        U_init = U_new;
+    
+        % Keep track of each iteration
+        iterations = iterations + 1;
+    end
+    caseIterations(numCases+1) = iterations;
+    numCases = numCases + 1;
+
+    U_solution = U_new;
+
+    % Reset conditions
+    U_init = U_original;
+    U_new = U_init;
+    maxError = 1;
+    iterations = 0;
+end
+
+% Plot the solution
+figure
+contourf(U_solution)
+title("Laplace equation solution")
+xlabel("x")
+ylabel("y")
+colorbar

FiniteDifferenceMethods/exercise/.gitkeep

@@ -0,0 +1,66 @@
+% Domain setup
+L = 1;
+H = 1;
+
+% Points in the domain
+nx = 22;
+ny = 22;
+
+% Grid resolution
+dx = L/nx;
+dy = H/ny;
+
+% Initializing solution
+U_init = zeros(ny,nx);
+
+% Boundary conditions
+U_init(:,1) = 10;
+U_init(:,nx) = 5;
+U_init(1,:) = 4;
+U_init(ny,:) = 1;
+
+U_sol = U_init;
+
+% Build the coefficient matrix
+% Use MATLAB function spdiag() and feed in matrix coefficients
+n = nx-2;
+e = ones(n,1);
+T = ;
+
+% Kronecker Delta product
+% Use MATLAb function kron, and take the prudct of I(x)T + T(x)I
+A = ;
+
+% Build the right hand side
+% This matrix will contain your boundary conditions and is equal to
+% the negative value of your stencil.
+for j=2:n+1
+    for i=2:n+1
+        rhs(j-1,i-1) = ;
+    end
+end
+
+% Reshape the rhs to a 1D vector
+% Use MATLAB reshape() to [n*n, 1]
+
+
+% Solve the equation x = A \ rhs
+
+
+% Reshape the matrix back to a n by n grid
+
+
+% Add back in the boundary conditions
+for j=2:n+1
+    for i=2:n+1
+        
+    end
+end
+
+%% Plot the solution
+figure
+contourf(U_sol)
+title("Laplace equation solution")
+xlabel("x")
+ylabel("y")
+colorbar

FiniteDifferenceMethods/exercise/hyperbolic_ex1.m

@@ -0,0 +1,68 @@
+% Domain setup
+L = 2;
+H = 1;
+
+% Points in the domain
+nx = 44;
+ny = 22;
+
+% Grid resolution
+dx = L/nx;
+dy = H/ny;
+
+% Initializing solution
+U_init = zeros (ny,nx);
+
+% Boundary conditions
+U_init(:,1) = 10;
+U_init(:,nx) = 5;
+U_init(1,:) = 4;
+U_init(ny,:) = 1;
+
+U_sol = U_init;
+
+% Build the coefficient matrix
+n_x = nx-2;
+n_y = ny-2;
+
+ex = ones(n_x,1);
+ey = ones(n_y,1);
+
+% Need Tx and Ty individually now
+Tx = ;
+Ty = ;
+
+% Kronecker Delta product
+% Product between I_x(x)Ty + Tx(x)I_y
+A = ;
+
+% Build the right hand side
+for j=2:n_y+1
+    for i=2:n_x+1
+        rhs(j-1,i-1) = ;
+    end
+end
+
+% Reshape the rhs to a 1D vector 
+
+
+% Solve the equation x = A \ rhs
+x = A\rhs;
+
+% Reshape the matrix back to a n by n grid
+
+
+% Add back in the boundary conditions
+for j=2:n_y+1
+    for i=2:n_x+1
+        
+    end
+end
+
+%% Plot the solution
+figure
+contourf(U_sol)
+title("Laplace equation solution")
+xlabel("x")
+ylabel("y")
+colorbar

FiniteDifferenceMethods/exercise/laplace_ex1.m

@@ -0,0 +1,49 @@
+% Define the parameters
+L = 1;             % Length of the rod
+T = 1;             % Total time
+N = 100;           % Number of spatial grid points
+alpha = 0.02;      % Thermal diffusivity constant
+
+for k=1:5
+    % Refine spatial resolution
+    if k > 1
+        N = N*2;
+    end
+
+    dx = L/N;
+    dt = 0.5*dx*dx/alpha;
+
+    M = int32(T/dt);
+
+    x = linspace(0,L,N+1);
+
+    % Initialize the temperature matrix
+    u = zeros(N+1,M+1);
+
+    % Set the initial condition
+    u(:,1) = sin(pi*x);
+    
+    time = 0.0;
+
+    for j = 1:M
+        for i = 2:N
+            u(i,j+1) = u(i,j) + alpha*dt/dx^2 * (u(i+1,j) - 2*u(i,j) + u(i-1,j));
+        end
+    
+        time = time + dt;
+        exact = sin(pi*x)*exp(-alpha*(pi)^2*time);
+    
+        Error = exact' - u(:,j);
+        totalError(j) = sum(Error);
+        
+    end
+    
+    % Error
+    E(k) = norm(totalError,inf);
+end
+
+%%
+figure
+plot(E,'o-')
+xlabel("Delta x")
+ylabel("Error")