initial commit

demo.m

+close all; clear all; clc;
+%% PAREMETERS
+x_left = 0.0;
+x_right = 1.0;
+rho = 1.0;
+c_p = 1.0;
+k   = 1.0;
+T_0 = 0.0;
+
+dt = 0.01;
+num_els = 100;
+num_nodes = num_els + 1;
+num_time_steps = 5;
+
+wait_time_plots = 1;%seconds
+
+%% PAREMETERS
+Xs = linspace(x_left,x_right,num_nodes);
+
+%% INITIAL CONDITION
+f0  = @(x) sin(16*x).^2 + (0.5 +0.5*cos(9*x));
+u_np1 = f0(Xs');
+
+ylims = [min(u_np1), max(u_np1)];
+
+fig = figure;
+plot(Xs, u_np1);
+ylim(ylims);
+pause(wait_time_plots);
+hold on
+
+%% TIME LOOP
+for time_iter=1:num_time_steps
+  u_n = u_np1;
+  % ASSEMBLY
+  A = sparse(num_nodes,num_nodes);
+  M = sparse(num_nodes,num_nodes);
+  b = zeros([num_nodes,1]);
+  for ielem=1:num_els
+    inodes_global = [ielem, ielem+1];
+    xr = Xs(ielem+1);
+    xl = Xs(ielem);
+    h  = xr - xl;
+    M_loc = h / 6.0 * [2.0, 1.0; 1.0, 2.0];
+    M_loc = M_loc * rho * c_p;
+    A_loc = 1 / h * [1.0, -1.0; -1.0, 1.0];
+    A_loc = A_loc * k;
+    M(inodes_global,inodes_global) = M(inodes_global,inodes_global) + M_loc;
+    A(inodes_global,inodes_global) = A(inodes_global,inodes_global) + A_loc;
+  end
+  % SOLVE
+  b = b + (1 / dt) * M * u_n;
+  u_np1 = mldivide((M / dt + A), b);
+  clf(fig);
+  plot(Xs, u_np1);
+  ylim(ylims);
+  pause(wait_time_plots)
+end