1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
% Coursework 17 Q4
% Using dx and dt as h and k are confusing
% h = dx
% k = dt
dx = input('Input position step (for example 0.01): '); % play about with this to get resolution
% Calulate maximum dt to maintain stability, based on the tailor expansion.
dt = dx^2/2;
tfin = input('Input the time you wawnt to end the simulation (for example 0.1): ');
lines = input('How many lines across the time range would you like to plot (for example 10): ');
% v = dx/(dt^2); %redundant
% Create x and t for plotting in the array
x = 0:dx:1;
t = 1:dt:tfin+1;
% Initialization of temperature matrix advancing in time (rows) and space (columns)
u = zeros(length(t),length(x));
u(1,:) = 0;
u(1,length(x)) = 0;
% Initial condition
for i = 1:length(x)
if x(i) <= 0.5
u(1,i) = 2*x(i);;
else
u(1,i) = 2*(1-x(i));
end
end
for m = 1:length(t)
% Set boundaries
u(m+1,1) = 0.5 * m /length(t);
u(m+1,length(x)) = 0.5 * m /length(t);
for j = 2:(length(x)-1)
% multiply out (1-2v) and factorise out v
u(m+1,j) = u(m,j) + ((dt/(dx^2))*(u(m,j+1) - 2*u(m,j) + u(m,j-1)));
end
end
figure;
hold on;
for i = 1:round(length(t)/lines):length(t)
plot(x,u(i,:),'.');
end
|
