Question 7.2: Let’s return to Example 7.1 and use explicit Euler to integr...
Let’s return to Example 7.1 and use explicit Euler to integrate the system of equations with α = 1 on the domain x ∈ [0, 1].
If we want to solve Eq. (7.2.11) using explicit Euler, the formula for propagating the solution in time becomes
\frac{dT}{dt} = f(T), T(0) = T_0 (7.2.11)
T^{(k+1)} = T^{(k)} + hf(T^{(k)}) (7.2.30)
We first need to modify Program 4.1 to solve this problem and output a vector of values T_i at each desired time step:
1 function [xout,yout] = ivp_explicit_euler(x0,y0,h,n_out,i_out)
2 % x0 = initial value for x
3 % y0 = initial value for y
4 % h = step size
5 % i_out = frequency to output x and y to [xout, yout]
6 % n_out = number of output points
7 % requires the function getf(x,y)
8 xout = zeros(n_out+1,1); xout(1) = x0;
9 yout = zeros(n_out+1,length(y0)); yout(1,:) = y0;
10 x = x0; y = y0;
11 for j = 2:n_out+1
12 for i = 1:i_out
13 y=y+h*getf(y);
14 x = x + h;
15 end
16 xout(j) = x;
17 yout(j,:) = y;
18 end
The changes required for a system of variables are in lines 8 and 9, where we set up the output, and lines 16 and 17, where we put in each output vector. We also need to write the program for computing the function f(T^{(k)}):
1 function f = getf(y)
2 n = length(y);
3 dx = 1/(n-1);
4 f = zeros(n,1);
5 f(1) = 0;
6 for i = 2:n-1
7 f(i) = y(i+1) - 2*y(i) + y(i-1);
8 end
9 f(n) = 0;
10 f = f/dx^2;
Since we are not passing the number of nodes to the function, we determine the value of n in line 2 by looking at the number of entries in the vector y. We then compute the value of Δx on line 3. When we fill in the entries in f(T^{(k)}), notice that we do not multiply each time by α/ Δx^2. Rather, we do this multiplication once at the end of the function in line 10.
Figure 7.2 shows the trajectory of the solution for a small value of h and a larger value of h. For the small value in Fig. 7.2a, we see that the solution is stable and we get the expected path to the steady state. However, as we increase h in Fig. 7.2b, the system eventually becomes numerically unstable. For this particular value of h the instability is gradual. If we make h even larger, the solution blows up very quickly.
