On the 17th of April from 18:00 to 20:00, We will be upgrading GitLab. During this time the service will be unavailable, apologies for any inconvenience
"The code is taking 100 iterations without stopping. **Add a condition to the code above to stop the loop once the solution is good enough**, i.e., when the solution is closer than $\\epsilon = 1e-6$ to the exact solution. How many iterations does it take now to converge? \n",
"The code is taking 100 iterations without stopping. **Add a condition to the code above to stop the loop once the solution is good enough**, i.e., when the solution is closer than $\\epsilon = 10^{-6}$ to the exact solution. How many iterations does it take now to converge? \n",
"\n",
"</p>\n",
"</div>"
...
...
@@ -301,7 +301,7 @@
"metadata": {},
"outputs": [],
"source": [
"def g(y_iplus1,y_i, t_iplus1):\n",
"def g(y_iplus1,y_i, t_iplus1):\n",
" return YOUR_CODE_HERE\n",
"\n",
"def g_der(y_iplus1):\n",
...
...
@@ -309,7 +309,7 @@
"\n",
"\n",
"# Define parameters:\n",
"dt = .3\n",
"dt = .25\n",
"t_end = 10\n",
"t = np.arange(0,t_end+dt,dt)\n",
"\n",
...
...
@@ -331,7 +331,7 @@
" y_IE[i+1] = YOUR_CODE_HERE # Initial guess\n",
" for j in range(200):\n",
" y_IE[i+1] = YOUR_CODE_HERE\n",
" if np.abs(g(y_IE[i+1],y_IE[i],t[i+1])) < 1e-6:\n",
" if np.abs(g(y_IE[i+1],y_IE[i],t[i+1])) < 1e-6:\n",
" break\n",
" \n",
" if j >= 199:\n",
...
...
@@ -529,17 +529,17 @@
"source": [
"# Initial conditions\n",
"\n",
"T_left = YOUR_CODE_HERE\n",
"T_right = YOUR_CODE_HERE\n",
"T_initial = YOUR_CODE_HERE\n",
"T_left = YOUR_CODE_HERE # Temperature at left boundary\n",
"T_right = YOUR_CODE_HERE # Temperature at right boundary\n",
This assignment contains two parts: treating non-linear ODEs and treating the diffusion equation (PDE).
%% Cell type:code id:453992c1 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import interact
```
%% Cell type:markdown id:0933143e tags:
## Part 1: Solving Non-linear ODEs
In task 1 you will solve first a very simple equation using Newton-Rhapson to understand exactly how to implement it. Task 2 treats the solution of a non-linear ODE w.r.t. time, first with Explicit Euler and then with Implicit Euler. The latter will require again Newton-Rhapson to find the solution.
Implement your equations $g(x)$ and $g'(x)$ in the code below, as well as the Newton-Rhapson expression inside the loop. Test the code with the initial guess of $x=10$.
</p>
</div>
%% Cell type:code id:25bcec4e tags:
``` python
import numpy as np
def g(x):
return YOUR_CODE_HERE
def g_der(x):
return YOUR_CODE_HERE
x = YOUR_CODE_HERE
for j in range(100):
x = YOUR_CODE_HERE
# Next task will go here
print("The solution found is ", x, " it took " ,j , " iterations to converge.")
The code is taking 100 iterations without stopping. **Add a condition to the code above to stop the loop once the solution is good enough**, i.e., when the solution is closer than $\epsilon = 1e-6$ to the exact solution. How many iterations does it take now to converge?
The code is taking 100 iterations without stopping. **Add a condition to the code above to stop the loop once the solution is good enough**, i.e., when the solution is closer than $\epsilon = 10^{-6}$ to the exact solution. How many iterations does it take now to converge?
Change the intial guess to $x=0.01$, which is closer to the exact solution than the initial guess in the previous task. How many iterations does it take to converge? Explain the difference.
Implement the Explicit Euler and Implicit Euler schemes by filling the lines of code below:
- a) Code the functions g() and g'()
- b) Implement Explicit Euler
- c) Implement Implicit Euler. Tip: use as initial guess the value of the previous time step.
- d) Use a dt = 0.25s
</p>
</div>
%% Cell type:code id:57c0662a tags:
``` python
def g(y_iplus1,y_i, t_iplus1):
def g(y_iplus1,y_i, t_iplus1):
return YOUR_CODE_HERE
def g_der(y_iplus1):
return YOUR_CODE_HERE
# Define parameters:
dt = .3
dt = .25
t_end = 10
t = np.arange(0,t_end+dt,dt)
y_EE = np.zeros(t.shape)
y_IE = np.zeros(t.shape)
# Define Initial Conditions
y_EE[0] = YOUR_CODE_HERE
y_IE[0] = YOUR_CODE_HERE
# Perform time-integration
newtonFailed = 0
for i in range(0, len(t)-1):
# Forward Euler:
y_EE[i+1] = YOUR_CODE_HERE
# Backward Euler:
y_IE[i+1] = YOUR_CODE_HERE # Initial guess
for j in range(200):
y_IE[i+1] = YOUR_CODE_HERE
if np.abs(g(y_IE[i+1],y_IE[i],t[i+1])) < 1e-6:
if np.abs(g(y_IE[i+1],y_IE[i],t[i+1])) < 1e-6:
break
if j >= 199:
newtonFailed = 1
# Plotting the solution
plt.plot(t, y_EE, 'r', t, y_IE, 'g--')
if newtonFailed:
plt.title('Nonlinear ODE with dt = ' + str(dt) + ' \nImplicit Euler did not converge')
else:
plt.title('Nonlinear ODE with dt = ' + str(dt))
plt.xlabel('t')
plt.ylabel('y')
plt.gca().legend(('Explicit','Implicit'))
plt.grid()
plt.show()
```
%% Cell type:markdown id:6b6d9964 tags:
## Part 2: Diffusion Equation in 1D
The 1-D diffusion equation reads $$\frac{\partial u}{\partial t}=v\frac{\partial^2 u}{\partial x^2}$$
where $u$ is a continuous function in space and time, $v$ is a constant and often referred to as the **diffusivity coefficient**, giving rise to the name 'diffusion equation'. This is a Partial Differential Equation which independent variable, $u$, varies on space and time. This equation is virtually present across all fields of civil engineering and science. Here, we use it to represent the temperature on a rod (see the sketch below).
Unlike the problem of Wednesday, here there is no exchange of heat with the ambient and the temperature evolves in time. The temperature initially is uniform along the rod, equal to $7°C$. Then it is heated at both ends with the temperatures shown in the figure.
The problem is schematized as a one-dimensional $0.3 m$ steel rod of with a diffusivity coefficient of $4e-6 m^2/s$. Run the simulation for $10,000 s$ to see the progression of the temperature through the model. Start with $200$ time steps and use 15 points to represent the rod.
Solve the diffusion equation using Central Differences in space and Forward Differences in time. You will do this step by step in the following sub-tasks.
For convenience we write here the diffusion equation with the temperature variable:
Draw a grid of 6 points with subindexes. Although your actual grid will be much larger, 6 points are enough to visualize the procedure. The initial condition states that the temperature of the rod is $7^o$ C. Does that mean that one single point of your grid is initialized or all of them?
Now, the differential equation needs to be expressed in algebraic form using central differences in space AND forward differences in time. **Start by just transforming the PDE into a first-order ODE by ONLY applying Central Differences to the spatial derivative term.**
Before applying Forward Differences in time to the equation. You need to add a superscript to the notation that indicates the time step: $T^j_i$. So, $i$ indicates the spatial location and $j$ the time location. For example, $T^0_2$ indicates the temperature at the node $i=2$ and at the initial moment $j=0 (t=0)$.
***Hint***: To express in a general form a node of over the next time step, you can write $T^{j+1}_i$.
**Apply Forward Differences to the equation to obtain an algebraic expression.**
Define the initial conditions and the boundary conditions. **Fill in the missing parts of the code.**
We define a 2-dimensional Numpy array <code>T</code> where the first index, <code>j</code>, represents time and the second index, <code>i</code>, represents space, for example: <code>T[j, i]</code>. Initialize <code>T</code> with a matrix of zeros.
Write in paper the equations that come out from your algebraic representation of the diffusion equation, solving for the unknowns. Use it then to write the matrix A, the unknown vector T and vector b. As in the workshop and textbook, the <code>A</code> matrix consists only of the unknowns in the problem.
Use this code cell if you would like to verify your numerical implementation. For example, to visualize the temperature profile at different time steps.
where L refers to the last point of the rod. Put your whole code together in a single cell. Copy the code cells from task 3.5 until task 3.8. Modify the right boundary condition as stated above, the period is 6000 seconds.
Solve the diffusion equation using Central Differences in space but **now with Backward Differences in time**. You will do this step by step just as before.
Draw the stencils (two in total) of this equation when solving it with Central Differences in space and **Forward Differences in time** and when solving it with Central Differences in space and **Backward Differences in time**.
Now, the differential equation needs to be expressed in algebraic form using central differences in space and forward differences in time. **Start by just transforming the PDE into a first-order ODE by ONLY applying Central Differences to the spatial derivative term.**
Write in paper the equations that come out from your algebraic representation of the diffusion equation, solving for the unknowns. Use it then to write the matrix A, the unknown vector T and vector b.
Copy the code of task 4 and make sure to use the Dirichlet conditions of task 3. Implement the Implicit scheme by modifying the code of how the matrix A and vector b are built.
