Problem Set Week 4

Yada Pruksachatkun

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
libwebsockets       par2                thefuck             youtube-dl
==> Renamed Formulae
recipes -> gnome-recipes

Warning: ffmpeg-3.2.4 already installed

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; 
}

window.implicit_ADI_diff = implicit_ADI_diff;
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

video = io.open('/Users/yadapruksachatkun/864.17/people/yada.pruksachatkun/homework/video-1489867333.mp4', 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))
Out[137]: