# Problem Set Week 4¶

In [24]:
!pip3 install matplotlib
!brew install ffmpeg

Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/site-packages
Requirement already satisfied: pytz in /usr/local/lib/python3.6/site-packages (from matplotlib)
Requirement already satisfied: six>=1.10 in /usr/local/lib/python3.6/site-packages (from matplotlib)
Requirement already satisfied: pyparsing!=2.0.0,!=2.0.4,!=2.1.2,!=2.1.6,>=1.5.6 in /usr/local/lib/python3.6/site-packages (from matplotlib)
Requirement already satisfied: numpy>=1.7.1 in /usr/local/lib/python3.6/site-packages (from matplotlib)
Requirement already satisfied: python-dateutil in /usr/local/lib/python3.6/site-packages (from matplotlib)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/site-packages (from matplotlib)
Updating Homebrew...
==> Auto-updated Homebrew!
Updated 1 tap (homebrew/core).
==> Updated Formulae
azure-cli           mosquitto           sonarqube           ttyd
==> Renamed Formulae
recipes -> gnome-recipes



8.1)

In [119]:
from IPython.display import Image
Image(filename='files/Week4Problem11.jpg')

Out[119]:

8.2

In [120]:
Image(filename='files/Week4Problem12.jpg')

Out[120]:
In [2]:
%%javascript
var abs = Math.abs;

function array_fill(i, n, v) {
var a = [];
for (i= 0; i < n; i++) {
a.push(v);
}
return a;
}
window.array_fill = array_fill;
function gauss(A, x) {
//TODO: Can implement this in numpy, slice and dice
var i, k, j;

for (i=0; i < A.length; i++) {
A[i].push(x[i]);
}
var n = A.length;

for (i=0; i < n; i++) {
var maxEl = abs(A[i][i]),
maxRow = i;
for (k=i+1; k < n; k++) {
if (abs(A[k][i]) > maxEl) {
maxEl = abs(A[k][i]);
maxRow = k;
}
}

for (k=i; k < n+1; k++) {
var tmp = A[maxRow][k];
A[maxRow][k] = A[i][k];
A[i][k] = tmp;
}

for (k=i+1; k < n; k++) {
var c = -A[k][i]/A[i][i];
for (j=i; j < n+1; j++) {
if (i===j) {
A[k][j] = 0;
} else {
A[k][j] += c * A[i][j];
}
}
}

}
var x = array_fill(0, n, 0);
for (i = (n-1);  i > -1; i--) {
x[i] = A[i][n]/A[i][i];
for (k=i-1; k > -1; k--) {
A[k][n] -= A[k][i] * x[i];
}
}
return x;
}
window.gauss = gauss;
function test_two() {
var A = [[1, 1], [2, 1]];
var x = [10, 16];
var result = gauss(A, x);
console.log(result);

}

test_two();


a)

In [74]:
%%javascript

function implicit_diff(D, sol_x, dt, dx, k) {
var A = [];
var local = [];
local.push(1)
var i;
var j;
var a;
for(i= 0; i < 499; i++) {
local.push(0);
}
A.push(local);
var alpha = (D*dt)/Math.pow(dx,2);
var first = -1*alpha;
var sec = 1+(2*alpha);
var third = alpha;
for(i= 0; i < 498; i++) {
local = []
for (j= 0; j < i; j++) {
local.push(0);
}
local.push(first);
local.push(sec);
local.push(third);
for (j = i+3 ; j < 500; j++) {
local.push(0);
}
A.push(local);

}
local = []
for(j= 0; j < 499; j++) {
local.push(0);
}
local.push(1);
A.push(local);

for (a = 0; a < k; a++) {
console.log("at time")
console.log(a);
sol_x = window.gauss(A, sol_x);

console.log(sol_x);
}
return;
}

window.implicit_diff = implicit_diff;


The diffusion equation should be unconditionally stable.

In [40]:
%%javascript
var t = [1, 0.5, 0.1];
var i;
var sol_x = window.array_fill(0, 249, 0);
sol_x.push(1);
for (i = 250; i< 500; i++) {
sol_x.push(0);
}

for (i = 0; i < 3; i++){
var dt = t[i]
console.log("for this value of dt")
console.log(dt);
window.implicit_diff(1, sol_x, dt, 1, 2);
}
// Here, you can now graph this.


6.3)

In [ ]:
import math
import numpy as np
# N is num iterations
def explicit_diff(D, N, iter_time, dx, dt):
F = dt/math.pow(dx,2)
u_n = ([0]*249) + [1] + ([0]*250)  # mesh points in space
u = [0]*500
final = []
final.append(u)
# Compute u at inner mesh points
for j in range(iter_time):
for i in range(1, (N-1)):
u[i] = u_n[i] + F*(u_n[i-1] - 2*u_n[i] + u_n[i+1])
print("At t=")
print(j)
print(u)
final.append(u)
u[0] = 0;
u[N-1] = 0
u_n = u
u = [0]*500
return final

dt = [1 , 0.5, 0.1]
for i in dt:
explicit_diff(1, 500, 10, 1, i)

# Now plot both these things.


I couldn't get the animations to work locally on Jupyter, so i put it on Youtube https://www.youtube.com/watch?v=k5mU1icwneA Here's the code for the animation

In [39]:
# ADI Method

In [128]:
%%javascript
// n by n matrix
function convert(A, alpha,dt, dx) {
// boundary conditions of 0
var i;
A[0] = 0;
A[A.length - 1] = 0;
for (i = 1; i < (A.length - 1); i++) {
A[i] = A[i] - (((alpha/2)*Math.pow((dt/dx), 2))*(A[i+1] -(2*A[i]) + A[i-1]));
}
return A;
}
// TO DO: Clena up, sinc ethis is similar to implicit_diff
function implicit_ADI_diff(D, sol_x, dt, dx, k) {
// console.log("here")
var A = [];
var local = [];
local.push(1)
var j;
var a;
var i;
for(i= 0; i < 399; i++) {
local.push(0);
}
A.push(local);
var alpha = (D*dt)/Math.pow(dx,2);
var first = -1*alpha;
var sec = 1+(2*alpha);
var third = alpha;
var i;
for(i= 0; i < 398; i++) {
local = []
for (j= 0; j < i; j++) {
local.push(0);
}
local.push(first);
local.push(sec);
local.push(third);
for (j = i+3 ; j < 400; j++) {
local.push(0);
}
A.push(local);
}
local = []
for(j= 0; j < 399; j++) {
local.push(0);
}

local.push(1);
A.push(local);
var final = []

for (a = 0; a < k; a++) {
//  console.log("now predicting stpes");
var half_step_x = window.gauss(A, sol_x);
half_step_x = convert(half_step_x, alpha, dt, dx)
var next_step = window.gauss(A, half_step_x);
final.push(next_step); //here, push a nxn next step
//  console.log("now predicting stpes");
}
return final;
}


In [129]:
%%javascript

var j;
var sol_x = []
for (j = 0; j< 160000; j++) {
var ran = Math.random();
sol_x.push(ran);
}
var z;
var t = [1, 0.5, 0.1];
var dt = t[0]
console.log("for this value of dt")
console.log(dt);
var animation = window.implicit_ADI_diff(1, sol_x, dt, 1, 10);
var json_string = JSON.stringify(animation);
var kernel = IPython.notebook.kernel;
console.log(json_string)

In [137]:
#Having a little trouble graphing 3D animations, can't mport matplotlib axes3D
# So instead here's a video of one runs, with 10 timeframes, and a 20x20 grid (400 sites) and dt=1

import io
import base64
from IPython.display import HTML