-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathApproxSolutionImplitEuler.m
111 lines (97 loc) · 2.89 KB
/
ApproxSolutionImplitEuler.m
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
function [z] = ApproxSolutionImplitEuler(Nx, Lx, Nt, ht, Asp, Bsp, input, z0)
%ApproxSolutionImplitEuler Computes and plot the solution to the full space
%semi-discretized problem, using the implicit euler scheme.
font = 15;
% Space axis
xaxis = linspace(0, Lx, Nx+2);
%int_mesh_xaxis=meshgrid(x1axis(1, 2:Nx+1), x1axis(1, 2:Nx+1));
[x1_axis, x2_axis] = meshgrid(xaxis, xaxis);
% Initial condition: heat the middle of the square
z = zeros(Nx^2, Nt+2); % rows<->space columns<->time
z(:,1) = reshape(z0', Nx^2, 1);
% Matrix to invert
M = sparse(1:Nx^2, 1:Nx^2, 1/ht) - Asp;
[L, U, P, Q] = lu(M);
% Implict Euler Scheme
for k=1:Nt+1
fk = 1/ht*z(:, k) + Bsp*input(:, k+1);
z(:, k+1) = Q*(U\(L\(P*fk))); % z(:, k+1) = M\fk;
end
%test
% zfull = zeros(Nx+2, Nx+2);
% test = z(:, 1); test(1, 1)=0.03; test(Nx, 1)=0.04; test(Nx^2-Nx+1, 1)=0.05; test(Nx^2, 1)=0.06;
% zreshape=reshape(test, Nx, Nx)';
% zfull(2:Nx+1, 2:Nx+1)=zreshape;
% surf(x1_axis, x2_axis, zfull);
% zlim([0 0.06]);
% xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
% title('Approximate solution of the PDE, at fixed time t')
% for k=2:Nt+2
% test = z(:, k); test(1, 1)=0.03; test(Nx, 1)=0.04; test(Nx^2-Nx+1, 1)=0.05; test(Nx^2, 1)=0.06;
% zreshape=reshape(test, Nx, Nx)';
% zfull(2:Nx+1, 2:Nx+1)=zreshape;
% surf(zfull);
% zlim([0 0.06]);
% drawnow;
% pause(0.1);
% end
%%% To plot :
figure;
zfull = zeros(Nx+2, Nx+2);
zreshape=reshape(z(:, 1), Nx, Nx)';
zfull(2:Nx+1, 2:Nx+1) = zreshape;
surf(x1_axis, x2_axis, zfull);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, at time t')
liste = [5 15 25 35];
save1 = zeros(Nx+2, Nx+2); save2 = zeros(Nx+2, Nx+2);
save3 = zeros(Nx+2, Nx+2); save4 = zeros(Nx+2, Nx+2);
for k=2:Nt+2
zreshape=reshape(z(:, k), Nx, Nx)';
zfull(2:Nx+1, 2:Nx+1) = zreshape;
if k==liste(1)
save1 = zfull;
end
if k==liste(2)
save2 = zfull;
end
if k==liste(3)
save3 = zfull;
end
if k==liste(4)
save4 = zfull;
end
surf(x1_axis, x2_axis, zfull);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, at time t')
set(gca,'fontsize',font)
drawnow;
pause(0.1);
end
figure;
surf(x1_axis, x2_axis, save1);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, iter=5')
set(gca,'fontsize',font)
figure;
surf(x1_axis, x2_axis, save2);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, iter=15')
set(gca,'fontsize',font)
figure;
surf(x1_axis, x2_axis, save3);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, iter=25')
set(gca,'fontsize',font)
figure;
surf(x1_axis, x2_axis, save4);
zlim([0 1]);
xlabel('x1'); ylabel('x2'); zlabel('z(t, x1, x2)')
title('Approximate solution, iter=35')
set(gca,'fontsize',font)
end