CS 131: Code for MP 3

Here is the code for MP 3

 

%%%heat equation
clear all
close all

x = 80; % size of the rod/beam
M = 500;
X = linspace(0,x,M);
t = 5; % time
c = 0.35; % heat equation constant

dx = X(2)-X(1);
dt = 0.01;
T = 0:dt:t;
T = logspace(0,t,length(T));
const = 0.35;
M = length(X) % number of mesh points
N = length(T) % number of time steps
r = c*dt/(dx**2)
% initialize U
A = zeros(1,M);
I = zeros(1,M);
A(M/2) = 10;
I(M/2) = 10;
for i = -10:10
A(M/2-i) = 100;
I(M/2-i) = 100;
endfor
% set matrix
d = zeros(1,M)’+(1-2*r);
upper = zeros(1,M-1)’+r;
Q = diag(d) + diag(upper,1) + diag(upper, -1);
U = zeros(N,M);
for i = 1:N
f = A.^2./(I + 3) – const*A;
A = A*Q – f;
g = A.^2 – I;
I = I*Q*const + g;
U(i,:) = A;
plot(X,A)
pause(0.00001)
endfor
pause
surf(X,T,U);
pause
imagesc(X,T,U);