Overview

This assignment is aimed to develop an understanding of the Ordinary Differential Equation (ODE). There will be two sections about cooling and heating scenerios, corresponding to the first-order and the second-order ODEs. Please go through the text that follows and perform all steps outlined therein.

Part 1: First-order ODE

In the study of heat transfer, Newton's law of cooling is a physical law which states that the rate of heat loss of a body is directly proportional to the difference in the temperatures between the body and its environment. It can be expressed in the form of ODE, as below:

$$\frac{dT}{dt}=-k(T - T_s)$$

where $T$ is the temperature of the object at time $t$, $T_s$ is the temperature of the surrounding and assumed to be constant, and $k$ is the constant that characterizes the ability of the object to exchange the heat energy (unit 1/s), which depends on the specific material properties.

Now, Let's consider an object with the initial temperature of 50°C in a surrounding environment with constant temperature at 20°C. The constant of heat exchange between the object and the environment is 0.5 $sec^{-1}$.

Overview

This assignment is aimed to develop an understanding of the Ordinary Differential Equation (ODE). There will be two sections about cooling and heating scenerios, corresponding to the first-order and the second-order ODEs. Please go through the text that follows and perform all steps outlined therein.

Part 1: First-order ODE

In the study of heat transfer, Newton's law of cooling is a physical law which states that the rate of heat loss of a body is directly proportional to the difference in the temperatures between the body and its environment. It can be expressed in the form of ODE, as below:

$$\frac{dT}{dt}=-k(T - T_s)$$

where $T$ is the temperature of the object at time $t$, $T_s$ is the temperature of the surrounding and assumed to be constant, and $k$ is the constant that characterizes the ability of the object to exchange the heat energy (unit 1/s), which depends on the specific material properties.

Now, Let's consider an object with the initial temperature of 50°C in a surrounding environment with constant temperature at 20°C. The constant of heat exchange between the object and the environment is 0.5 $s^{-1}$.

Task 1.1:

Suppose the considered period of time is long enough (bring it to steady state), what will be the final temperature of the object?

Write your answer in the following markdown cell.

20°C

This is because the body's temperature, $T$, will tend to the ambience temperature in the long term. The decrease will be abrupt at the beginning and will slow down as time passes $T(t=large) = T_s = 20°C$.

Next, let's evaluate the temperature of the object by checking it at a series of time points.

Task 1.2:

Write the algebraic representation of the ODE using Explicit Euler.

Answer for Task 1.2:

$$ T_{i+1} = T_{i}-k(T_{i} - T_s)dt$$

How to get there:

$$ \int_{t}^{t+1} \frac{dT}{dt} dt=\int_{t}^{t+1}-k(T - T_s)dt $$$$ \int_{t}^{t+1} T' dt=\int_{t}^{t+1}f(t)dt $$

FTC (fundamental theorem of calculus) for left side and using left rieman to approximate right term

$$ T_{t+1} - T_t = \Delta t f(t_i) $$$$ T_{t+1}= T_t + \Delta t f(t_i) $$$$ T_{t+1}= T_t - k(T_t - T_s) $$

or from forward finite difference

Task 1.3:

Compute the temperature evolution in the next 60 seconds.

Please complete the missing parts of the code in each step below, which is divided into 5 substeps (a through e).

Task 1.3a:

The time step of 5 seconds is constant. Discretize the time points, the solution vector $T$ and define the initial condition.

import numpy as np
import matplotlib.pyplot as plt

dt =5                          
t_end = 60  
Ts =  20       # [C] 
k = 0.5        # [s^-1]

t = np.arange(0,t_end+dt,dt)
n = len(t)
T = np.empty(n)
T[0] = 50

Task 1.3b:

Implement your time discretization and find the solution from $t=0$ until $t=60$ sec.

for i in range(n-1):
    T[i+1] = T[i]-k*(T[i]-Ts)*dt 
    
plt.plot(t, T, 'o-')
plt.xlabel('t (s)')
plt.ylabel('T (deg)')
plt.grid()

Task 1.3c:

Try different time steps to check the stability of the calculation. At which value the solution is stable?

Answer:

4 seconds is the limit between an unstable and a stable solution.

Task 1.3d:

Obtain the mathematical expression that proves your stability criteria.

Answer

$$ \frac{dT}{dt}= -k(T-T_s) $$$$ \int \frac{1}{(T-T_s)}dT=\int -k dt $$$$ ln|T-T_s| = -kt +C_1 $$$$ T-T_s = e^{-kt+C_1} $$

Since $C_1$ is an arbitrary constant we simply define a new constant $C_2$:

$$ T-T_s = C_2e^{-kt} $$$$ T= T_s +C_2e^{-kt} $$

This is an Explicit Euler scheme, thus is conditionally stable, with the criterion: $ |1-k*dt|<1$.

Solving for n time steps, the solution can be generalized as $T_n=T_0(1-k\Delta t)^n + k \Delta t T_s \left[ (1-k \Delta t)^{n-1} + (1-k \Delta t)^{n-2} + ... + (1-k \Delta t)^{0} \right]$

The criterion of stability is $ |1-k*dt|<1$. Otherwise, the solution would tend to grow instead of decaying. Then $|k*dt| < 2$, thus $dt < 4$.

Task 1.3e:

Now, discretize the equation using Implicit (Backward) Euler. Can you find a time step that makes the problem unstable?

dt = 10                          
t_end = 60  
Ts =  20       # [C] 
k = 0.5        # [s^-1]

t = np.arange(0,t_end+dt,dt)
n = len(t)
T = np.empty(n)
T[0] = 50     

for i in range(n-1):
    T[i+1] = (T[i]+k*Ts*dt)/(1+k*dt) 
    
plt.plot(t, T, 'o-')
plt.xlabel('t (s)')
plt.ylabel('T (deg)')
plt.grid()

Answer:

No, you cannot! This method is unconditionally stable. (a good thing to remember for the exam)

Part 2: Second-order ODE

The following 1D equation describes the steady state solution of the temperature along a pin that sticks out of a furnace. The rest of the pin is exposed to the ambient.

$$ \frac{d^2T}{dx^2} -\alpha(T-T_s)=0 $$

The ambient temperature is $T_s= 30^o$ C and the temperature at the wall is $250^o$ C. The length of the pin is 0.1m. Your grid should have a spatial step of 0.02 m. Finally, $\alpha=500$.

The solution includes the steps:

1. Use the Taylor series to obtain an approximation for the derivatives;
2. Discretize the equation;
3. Define parameters and grid;
4. Provide boundary conditions;
5. Build matrix with solution $AT=b$
6. Solve the matrix

Task 2.1:

This task has three parts: a) discretize the analytic expression into a system of equations using central differences, b) write the system of equations by hand (e.g., the matrices), and c) implement the discretized system of equations in a code cell.

Parts 2.1b and 2.1c do not need to be completed in order; in fact, it may be useful to go back and forth between the two in order to understand the problem.

Task 2.1a:

Discretize the analytic expression into a system of equations for a grid with 6 points using central differences.

Write your answer by hand using paper/tablet and include the image below.

Answer for 2.1a and 2.1b

$$ \frac{1}{\Delta x^2}(T_{i-1} - 2T_{i} + T_{i+1}) - \alpha T_{i} + \alpha T_s = 0 $$

For $i = 1$: The finite difference approximation for the first internal node is:

$$ \frac{1}{\Delta x^2}(T_0 - 2T_1 + T_2) - \alpha T_1 + \alpha T_s = 0 $$

Next, we multiply by $\Delta x^2$ and substitute the boundary condition for $T_0$. Notice that the only unknown solutions are $T_1$ and $T_2$ Rearranging the terms gives:

$$ -(2 + \alpha \Delta x^2) T_1 + T_2 = -\alpha T_s \Delta x^2 - T_0 $$

For $i = 2$: The finite difference approximation for the second internal node is:

$$ \frac{1}{\Delta x^2}(T_1 - 2T_2 + T_3) - \alpha T_2 + \alpha T_s = 0 $$

Rearranging the terms:

$$ T_1 - (2 + \alpha \Delta x^2) T_2 + T_3 = -\alpha T_s \Delta x^2 $$

For $i = 3$: The finite difference approximation for the third internal node is:

$$ \frac{1}{\Delta x^2}(T_2 - 2T_3 + T_4) - \alpha T_3 + \alpha T_s = 0 $$

Notice that now we can plug in our other bc $T_5$ and rearrange the terms:

$$ (T_3 - (2 + \alpha \Delta x^2) T_4 ) = -\alpha Ts \Delta x^2 -T_5 $$

Task 2.1b:

Write the system of equations by hand for a grid with 6 points. In other words, construct the matrix A and vectors T and b.

Write your answer by hand using paper/tablet and include the image below.

The discretized equations are presented in Ax=y form:

$$ \begin{bmatrix} -(2 + \alpha \Delta x^2) & 1 & 0 & 0 \\ 1 & -(2 + \alpha \Delta x^2) & 1 & 0 \\ 0 & 1 & -(2 + \alpha \Delta x^2) & 1 \\ 0 & 0 & 1 & -(2 + \alpha \Delta x^2) \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \\ T_3 \\ T_4 \end{bmatrix} = \begin{bmatrix} -\alpha Ts \Delta x^2 -T_0\\ -\alpha Ts \Delta x^2 \\ -\alpha Ts \Delta x^2\\ -\alpha Ts \Delta x^2 -T_5 \end{bmatrix} $$

Task 2.1c:

Implement the discretized system of equations in a code cell.

We have already done this for you! Your task is to read the cell and make sure you understand how the matrices are implemented. Reading the code should help you formulate the matrices in Task 2.1b.

This is very similar to the code in the book, except that the matrix size here adapts to the number of grid points that are implied by the dx value selected.

import numpy as np 
import matplotlib.pyplot as plt

Ts = 30
alpha = 500
dx=0.02

# grid creation
x = np.arange(0,0.1+dx,dx)
T = np.zeros(x.shape)
n=len(x)

# boundary conditions
T[0] = 250
T[-1] = Ts

# Building matrix A
matrix_element = -(2+dx**2*alpha)
A = np.zeros((len(x)-2,len(x)-2))
np.fill_diagonal(A, matrix_element)
A[np.arange(n-3), np.arange(1, n-2)] = 1  # Upper diagonal
A[np.arange(1, n-2), np.arange(n-3)] = 1  # Lower diagonal

# Building vector b
b_element = -dx**2*alpha*Ts
b = np.zeros(len(x)-2) + b_element
b[0] = b[0] - T[0]
b[-1] = b[-1] - T[-1]

# Solving the system

T[1:-1] = np.linalg.solve(A,b)

plt.plot(x,T,'*',label='Estimated solution')
plt.xlabel('x')
plt.ylabel('T')
plt.title('Estimated solution using Central Difference method')
plt.legend()
plt.show()

print(f'The estimated temperature at the nodes are: {[f"{temp:.2f}" for temp in T]} [C]')

The estimated temperature at the nodes are: ['250.00', '168.77', '115.29', '78.86', '52.21', '30.00'] [C]

Task 2.2:

This task will adapt the problem from 2.1 to incorporate one Dirichlet condition (at the left boundary) and one Neumann boundary condition (at the right boundary). We will accomplish this in three steps: a) writing the new matrix by hand, b) adapting the code from 2.1c, c) reflecting on what this represents physically.

Task 2.2a:

Write the system of equations by hand for a grid with 6 points, incorporating the Neumann condition at the right boundary.

Approximate the Neuman boundary $\frac{dT}{dx}=0$ by using the Backward difference for first order differential equation of first order accuracy.

Write your answer by hand using paper/tablet and include the image below.

The first elements for i = 1,2,3 are the same as above in task 5

The neuman bc is on the last node

$$ \frac{dT}{dx} = \frac{T_i - T_{i-1}}{\Delta x} $$

Neuman bc $$\frac{dT_5}{dx}=0 \approx \frac{T_5 - T_{4}}{\Delta x} $$ Therefore $T_5=T_4$

The only equation that changes is for $i = 4$:

$$ \frac{1}{\Delta x^2}(T_3-2T_4+T_5) - \alpha T_4 + \alpha Ts =0 $$

Plug in BC $T_5=T_4$

$$ \frac{1}{\Delta x^2}(T_3-2T_4+T_4) - \alpha T_4 + \alpha Ts =0 $$$$ T_3-(1 + \alpha \Delta x^2)T_4 = \alpha Ts $$

We can now move the equations in the Ax=y form:

$$ \begin{bmatrix} -(2 + \alpha \Delta x^2) & 1 & 0 & 0 \\ 1 & -(2 + \alpha \Delta x^2) & 1 & 0 \\ 0 & 1 & -(2 + \alpha \Delta x^2) & 1 \\ 0 & 0 & 1 & -(1 + \alpha \Delta x^2) \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \\ T_3 \\ T_4 \end{bmatrix} = \begin{bmatrix} -\alpha Ts \Delta x^2 -T_0\\ -\alpha Ts \Delta x^2 \\ -\alpha Ts \Delta x^2\\ -\alpha Ts \Delta x^2 \end{bmatrix} $$

Task 2.2b:

Now adapt the code from Task 2.1c and revise it to incorporate the Neumann boundary condition.

Copy and past the code from 2.1c below, then modify it.

Compared to Task 2.1c there are three lines of code that change.

import numpy as np
import matplotlib.pyplot as plt
 
Ts = 30
alpha = 500
dx=0.02
 
# grid creation
x = np.arange(0,0.1+dx,dx)
T = np.zeros(x.shape)
n=len(x)
 
# boundary conditions
T[0] = 250
 
 
# Building matrix A
matrix_element = -(2+dx**2*alpha)
A = np.zeros((len(x)-2,len(x)-2))
np.fill_diagonal(A, matrix_element)
A[np.arange(n-3), np.arange(1, n-2)] = 1  # Upper diagonal
A[np.arange(1, n-2), np.arange(n-3)] = 1  # Lower diagonal
 
A[-1,-1] = -(1+dx**2*alpha)  # the lower right corner of the matrix changes
 
# Building vector b
b_element = -dx**2*alpha*Ts
b = np.zeros(len(x)-2) + b_element
b[0] = b[0] - T[0]
b[-1] = b[-1]               # the vector b also changes
# the line above does not change the value of b[-1];
# it is kept in place to help compare to Task 2.1c
 
# Solving the system

T[1:-1] = np.linalg.solve(A,b)
T[-1] = T[-2]  # this line has been added
 
plt.plot(x,T,'*',label='Estimated solution')
plt.xlabel('x')
plt.ylabel('T')
plt.title('Estimated solution using Central Difference method with Neumann BC')
plt.legend()
plt.show()
 
print(f'The estimated temperature at the nodes are: {[f"{temp:.2f}" for temp in T]} [C]')

The estimated temperature at the nodes are: ['250.00', '174.84', '128.64', '102.18', '90.15', '90.15'] [C]

Task 2.2c:

Reflect on the difference between the problem solved in Task 2.1 in comparison to 2.2. How are we changing the physics of the problem being solved by changing the boundary condition? What does this mean in reality for the temperature distribution in the bar over time?

We have changed the boundaries from a temperature to the derivative of the temperature (with respect to space). This is like changing it from a constant temperature at the edges to a of heat going in or out of the system at one of the edges. In this case since the gradient is 0 it represents a perfectly insulated system. Think of it like a really good portable cooler keeping your beverages cold at the beach on a hot day.

Task 2.3:

Just as we did in Task 2.1, this task has three parts: a) discretize the analytic expression into a system of equations using forward differences, b) write the system of equations by hand (e.g., the matrices), and c) implement the discretized system of equations in a code cell.

Here we focus on Dirichlet conditions again.

Task 2.3a:

Discretize the analytic expression into a system of equations for a grid with 6 points using forward differences.

Write your answer by hand using paper/tablet and include the image below.

Forward difference for second order differential equation: $$ \frac{T_i-2T_{i+1}+T_{i+2}}{\Delta x^2}-\alpha(T_i-T_s) = 0 $$

For $i = 0$: The finite difference approximation for the first internal node is:

$$ \frac{1}{\Delta x^2}(T_0 - 2T_1 + T_2) - \alpha T_0 + \alpha T_s = 0 $$$$ -2T_1 +T_2 = -(1-\alpha \Delta x^2)T_0 - \alpha T_s \Delta x^2 $$

For $i = 1$: The finite difference approximation for the second internal node is:

$$ \frac{1}{\Delta x^2}(T_1 - 2T_2 + T_3) - \alpha T_1 + \alpha T_s = 0 $$

Rearranging the terms:

$$ (1 - \alpha \Delta x^2)T_1 - 2T_2 + T_3 = -\alpha T_s \Delta x^2 $$

For $i = 2$: The finite difference approximation for the second internal node is:

$$ \frac{1}{\Delta x^2}(T_2 - 2T_3 + T_4) - \alpha T_2 + \alpha T_s = 0 $$

Rearranging the terms:

$$ (1 - \alpha \Delta x^2)T_2 - T_3 + T_4 = -\alpha T_s \Delta x^2 $$

For $i = 3$: The finite difference approximation for the second internal node is:

$$ \frac{1}{\Delta x^2}(T_3 - 2T_4 + T_5) - \alpha T_3 + \alpha T_s = 0 $$

Rearranging the terms:

$$ (1 - \alpha \Delta x^2)T_3 - 2T_4 = -\alpha T_s \Delta x^2 -T_5 $$

Task 2.3b:

Write the system of equations by hand for a grid with 6 points. In other words, construct the matrix A and vectors T and b.

Write your answer by hand using paper/tablet and include the image below.

We can now move the equations in the Ax=y form:

$$ \begin{bmatrix} -2 & 1 & 0 & 0 \\ (1- \alpha \Delta x^2) & -2 & 1 & 0 \\ 0 & (1- \alpha \Delta x^2) & -2 & 1 \\ 0 & 0 & (1- \alpha \Delta x^2) & -2 \end{bmatrix} \begin{bmatrix} T_1 \\ T_2 \\ T_3 \\ T_4 \end{bmatrix} = \begin{bmatrix} -\alpha Ts \Delta x^2 - (1-\alpha \Delta x^2)T_0\\ -\alpha Ts \Delta x^2 \\ -\alpha Ts \Delta x^2\\ -\alpha Ts \Delta x^2 -T_5 \end{bmatrix} $$

Task 2.3c:

Implement the discretized system of equations in a code cell.

This time we did not do it for you! Copy the code from Task 2.1c and revise it to solve the system of equations using Forward Differences. Keep the Dirichlet conditions.

import numpy as np 
import matplotlib.pyplot as plt

Ts = 30
alpha = 500
dx=0.0004 # The grid size is reduced to (0.02^2 = 0.0004)

# grid creation
x = np.arange(0,0.1+dx,dx)
T = np.zeros(x.shape)
n=len(x)

# boundary conditions
T[0] = 250
T[-1] = Ts

# Building matrix A
matrix_element = 1-alpha*dx**2    # The matrix element changes (from CD)
A = np.zeros((len(x)-2,len(x)-2))
np.fill_diagonal(A, -2)            # different from CD     
A[np.arange(n-3), np.arange(1, n-2)] = 1  # Upper diagonal
A[np.arange(1, n-2), np.arange(n-3)] = matrix_element  # Lower d.: different from CD
print(A)

# Building vector b
b_element = -dx**2*alpha*Ts
b = np.zeros(len(x)-2) + b_element
b[0] = b[0] - matrix_element * T[0]
b[-1] = b[-1] - T[-1]

# Solving the system

T[1:-1] = np.linalg.solve(A,b)

plt.plot(x,T,'*',label='Estimated solution')
plt.xlabel('x')
plt.ylabel('T')
plt.title('Estimated solution using Forward Difference method')
plt.legend()
plt.show()

print(f'The x values =============================: {[f"{x_loc:.4f}" for x_loc in x]} [m]')
print(f'The estimated temperature at the nodes are: {[f"{temp:.2f}" for temp in T]} [C]')

ind = np.argmin(abs(x-0.02))
print(f'The temperature at x=0.02 is: {T[ind]:.2f} [C]')

[[-2.       1.       0.      ...  0.       0.       0.     ]
 [ 0.99992 -2.       1.      ...  0.       0.       0.     ]
 [ 0.       0.99992 -2.      ...  0.       0.       0.     ]
 ...
 [ 0.       0.       0.      ... -2.       1.       0.     ]
 [ 0.       0.       0.      ...  0.99992 -2.       1.     ]
 [ 0.       0.       0.      ...  0.       0.99992 -2.     ]]

The x values =============================: ['0.0000', '0.0004', '0.0008', '0.0012', '0.0016', '0.0020', '0.0024', '0.0028', '0.0032', '0.0036', '0.0040', '0.0044', '0.0048', '0.0052', '0.0056', '0.0060', '0.0064', '0.0068', '0.0072', '0.0076', '0.0080', '0.0084', '0.0088', '0.0092', '0.0096', '0.0100', '0.0104', '0.0108', '0.0112', '0.0116', '0.0120', '0.0124', '0.0128', '0.0132', '0.0136', '0.0140', '0.0144', '0.0148', '0.0152', '0.0156', '0.0160', '0.0164', '0.0168', '0.0172', '0.0176', '0.0180', '0.0184', '0.0188', '0.0192', '0.0196', '0.0200', '0.0204', '0.0208', '0.0212', '0.0216', '0.0220', '0.0224', '0.0228', '0.0232', '0.0236', '0.0240', '0.0244', '0.0248', '0.0252', '0.0256', '0.0260', '0.0264', '0.0268', '0.0272', '0.0276', '0.0280', '0.0284', '0.0288', '0.0292', '0.0296', '0.0300', '0.0304', '0.0308', '0.0312', '0.0316', '0.0320', '0.0324', '0.0328', '0.0332', '0.0336', '0.0340', '0.0344', '0.0348', '0.0352', '0.0356', '0.0360', '0.0364', '0.0368', '0.0372', '0.0376', '0.0380', '0.0384', '0.0388', '0.0392', '0.0396', '0.0400', '0.0404', '0.0408', '0.0412', '0.0416', '0.0420', '0.0424', '0.0428', '0.0432', '0.0436', '0.0440', '0.0444', '0.0448', '0.0452', '0.0456', '0.0460', '0.0464', '0.0468', '0.0472', '0.0476', '0.0480', '0.0484', '0.0488', '0.0492', '0.0496', '0.0500', '0.0504', '0.0508', '0.0512', '0.0516', '0.0520', '0.0524', '0.0528', '0.0532', '0.0536', '0.0540', '0.0544', '0.0548', '0.0552', '0.0556', '0.0560', '0.0564', '0.0568', '0.0572', '0.0576', '0.0580', '0.0584', '0.0588', '0.0592', '0.0596', '0.0600', '0.0604', '0.0608', '0.0612', '0.0616', '0.0620', '0.0624', '0.0628', '0.0632', '0.0636', '0.0640', '0.0644', '0.0648', '0.0652', '0.0656', '0.0660', '0.0664', '0.0668', '0.0672', '0.0676', '0.0680', '0.0684', '0.0688', '0.0692', '0.0696', '0.0700', '0.0704', '0.0708', '0.0712', '0.0716', '0.0720', '0.0724', '0.0728', '0.0732', '0.0736', '0.0740', '0.0744', '0.0748', '0.0752', '0.0756', '0.0760', '0.0764', '0.0768', '0.0772', '0.0776', '0.0780', '0.0784', '0.0788', '0.0792', '0.0796', '0.0800', '0.0804', '0.0808', '0.0812', '0.0816', '0.0820', '0.0824', '0.0828', '0.0832', '0.0836', '0.0840', '0.0844', '0.0848', '0.0852', '0.0856', '0.0860', '0.0864', '0.0868', '0.0872', '0.0876', '0.0880', '0.0884', '0.0888', '0.0892', '0.0896', '0.0900', '0.0904', '0.0908', '0.0912', '0.0916', '0.0920', '0.0924', '0.0928', '0.0932', '0.0936', '0.0940', '0.0944', '0.0948', '0.0952', '0.0956', '0.0960', '0.0964', '0.0968', '0.0972', '0.0976', '0.0980', '0.0984', '0.0988', '0.0992', '0.0996', '0.1000'] [m]
The estimated temperature at the nodes are: ['250.00', '247.99', '245.99', '244.01', '242.05', '240.11', '238.18', '236.27', '234.38', '232.50', '230.64', '228.80', '226.97', '225.16', '223.36', '221.58', '219.81', '218.06', '216.33', '214.61', '212.90', '211.21', '209.54', '207.88', '206.23', '204.60', '202.98', '201.37', '199.78', '198.21', '196.64', '195.09', '193.56', '192.03', '190.52', '189.03', '187.54', '186.07', '184.61', '183.16', '181.73', '180.31', '178.90', '177.50', '176.11', '174.74', '173.38', '172.03', '170.69', '169.36', '168.04', '166.74', '165.44', '164.16', '162.88', '161.62', '160.37', '159.13', '157.90', '156.68', '155.47', '154.27', '153.08', '151.90', '150.73', '149.57', '148.42', '147.28', '146.15', '145.02', '143.91', '142.81', '141.71', '140.63', '139.55', '138.48', '137.42', '136.37', '135.33', '134.30', '133.27', '132.26', '131.25', '130.25', '129.26', '128.27', '127.30', '126.33', '125.37', '124.42', '123.47', '122.54', '121.61', '120.68', '119.77', '118.86', '117.96', '117.07', '116.18', '115.30', '114.43', '113.56', '112.70', '111.85', '111.00', '110.16', '109.33', '108.50', '107.68', '106.87', '106.06', '105.26', '104.47', '103.68', '102.89', '102.12', '101.34', '100.58', '99.82', '99.06', '98.32', '97.57', '96.83', '96.10', '95.37', '94.65', '93.94', '93.22', '92.52', '91.82', '91.12', '90.43', '89.74', '89.06', '88.38', '87.71', '87.04', '86.38', '85.72', '85.07', '84.42', '83.78', '83.13', '82.50', '81.87', '81.24', '80.61', '80.00', '79.38', '78.77', '78.16', '77.56', '76.96', '76.36', '75.77', '75.18', '74.60', '74.01', '73.44', '72.86', '72.29', '71.72', '71.16', '70.60', '70.04', '69.49', '68.94', '68.39', '67.85', '67.30', '66.77', '66.23', '65.70', '65.17', '64.64', '64.12', '63.60', '63.08', '62.56', '62.05', '61.54', '61.03', '60.53', '60.02', '59.52', '59.02', '58.53', '58.03', '57.54', '57.05', '56.57', '56.08', '55.60', '55.12', '54.64', '54.17', '53.69', '53.22', '52.75', '52.28', '51.81', '51.35', '50.89', '50.42', '49.96', '49.51', '49.05', '48.59', '48.14', '47.69', '47.24', '46.79', '46.34', '45.90', '45.45', '45.01', '44.56', '44.12', '43.68', '43.24', '42.80', '42.37', '41.93', '41.50', '41.06', '40.63', '40.20', '39.77', '39.34', '38.91', '38.48', '38.05', '37.62', '37.19', '36.77', '36.34', '35.92', '35.49', '35.07', '34.64', '34.22', '33.80', '33.37', '32.95', '32.53', '32.11', '31.69', '31.26', '30.84', '30.42', '30.00'] [C]
The temperature at x=0.02 is: 168.04 [C]

Task 2.4:

How much finer does your grid has to be in the forward difference implementation to get a similar value at x = 0.02 as in the central difference implementation? Vary your dx.

You can test this above make the temperature the same at x=0.02

Bonus Task

The matrix inversion using numpy is one way to solve the system, another is the gauss_jordan method, written below, and another one is the sparse matrix-based method in the cell afterwards. Here, we will just have a brief comparison to see how these solvers perform when the matrix is large. Change dx to 0.0002 of the original code that solves the second degree ODE and test the time it takes by each method.

def gauss_jordan(A, b):
    """
    Solves the system of linear equations Ax = b using Gauss-Jordan elimination.
    
    Parameters:
    A (numpy.ndarray): Coefficient matrix (n x n).
    b (numpy.ndarray): Right-hand side vector (n).
    
    Returns:
    numpy.ndarray: Solution vector (x) if the system has a unique solution.
    """
    # Form the augmented matrix [A | b]
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)
    aug_matrix = np.hstack([A, b.reshape(-1, 1)])
    
    n = len(b)  # Number of rows (or variables)
    
    for i in range(n):
        # Partial pivoting to handle zero diagonal elements (optional, but more robust)
        max_row = np.argmax(np.abs(aug_matrix[i:, i])) + i
        if aug_matrix[max_row, i] == 0:
            raise ValueError("The matrix is singular and cannot be solved.")
        if max_row != i:
            aug_matrix[[i, max_row]] = aug_matrix[[max_row, i]]
        
        # Make the diagonal element 1
        aug_matrix[i] = aug_matrix[i] / aug_matrix[i, i]
        
        # Make all other elements in the current column 0
        for j in range(n):
            if j != i:
                aug_matrix[j] -= aug_matrix[j, i] * aug_matrix[i]
    
    # Extract the solution (last column of the augmented matrix)
    return aug_matrix[:, -1]

import time
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import spsolve

# Inverted matrix solution
start_time = time.time()
A_inv = np.linalg.inv(A)
T[1:-1] = A_inv @ b
time_used_0 = time.time() - start_time
print(f"The time used by direct matrix inversion solution is {time_used_0: 0.3e} sec")
assert np.allclose(np.dot(A, T[1:-1]), b), "Oops! The calculation is wrong.."


# Gauss-jordan solution
start_time = time.time()
u1 = gauss_jordan(A, b)
time_used_1 = time.time() - start_time
print(f"The time used by Gauss-jordan solution is {time_used_1: 0.3e} sec")
#Check if the solution is correct:
assert np.allclose(np.dot(A, u1), b), "Oops! The calculation is wrong.."

# Solution by a sparse matrix solver 
start_time = time.time()
A_sparse = csc_matrix(A)# Convert A to a compressed sparse column (CSC) matrix
u2 = spsolve(A_sparse, b)
time_used_2 = time.time() - start_time
print(f"The time used by the sparse matrix solver is {time_used_2: 0.3e} sec")
#Check if the solution is correct:
assert np.allclose(np.dot(A, u2), b), "Oops! The calculation is wrong.."

The time used by direct matrix inversion solution is  1.399e-03 sec
The time used by Gauss-jordan solution is  1.411e-01 sec
The time used by the sparse matrix solver is  1.286e-03 sec

Answer

The sparse matrix is fastest, but we see that the Numpy implementation is also very fast! problem still small; Gauss-Jordan slow because...???

End of notebook.

© Copyright 2023 MUDE Teaching Team TU Delft. This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.