summaryrefslogtreecommitdiff
path: root/coursework17/ex4.m
blob: b480258157dffca1787116cc3f0e169cd691fca3 (plain)
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
% Coursework 17 Q4
% Using dx and dt as h and k are confusing
% h = dx
% k = dt

clear

dx = 0.002
% 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 = 0.1
% tfin = input('Input the time you wawnt to end the simulation (for example 0.1): ');

lines = 10
% lines = input('How many lines across the time range would you like to plot (for example 10): ');

% 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 1
% 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

% % Initial condition 2
% for i = 1:length(x)
% %	u(1,i) = abs(sin(2*pi*x(i)));
% 	u(1,i) = sin(2*pi*x(i));
% end

% % Initial Condition 3 
%  for i = 0.25/dx:0.75/dx
%  	u(1,i) = 1;
%  end
%

 % Initial Condition 4
 for i = 1:0.25/dx
 	u(1,i) = -1;
 end
 for i = 0.25/dx+1:0.75/dx
 	u(1,i) = -1+(i*dx-0.25)*4;
 end
 for i = 0.75/dx+1:length(x)
 	u(1,i) = 1;
 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); 
    u(m+1,1) = 0;
    u(m+1,length(x)) = 0;
    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;
j = 0
for i = 1:round(length(t)/lines):length(t)
    j = j+1;
    plot(x,u(i,:),'.');
    legendInfo{j} = ['t = ' num2str(round(i*dt, 3))];
end
legend(legendInfo,'Location', 'southeast');
xlabel('Displacement');
ylabel('Temperature');