Working code

README.md

@@ -0,0 +1 @@
+Solve the boundary layer velocity profile using a finite difference method. Directly animate the result.

blflow.lyx

@@ -0,0 +1,165 @@
+#LyX 2.1 created this file. For more info see http://www.lyx.org/
+\lyxformat 474
+\begin_document
+\begin_header
+\textclass article
+\use_default_options true
+\maintain_unincluded_children false
+\language english
+\language_package default
+\inputencoding auto
+\fontencoding global
+\font_roman default
+\font_sans default
+\font_typewriter default
+\font_math auto
+\font_default_family default
+\use_non_tex_fonts false
+\font_sc false
+\font_osf false
+\font_sf_scale 100
+\font_tt_scale 100
+\graphics default
+\default_output_format default
+\output_sync 0
+\bibtex_command default
+\index_command default
+\paperfontsize default
+\use_hyperref false
+\papersize default
+\use_geometry false
+\use_package amsmath 1
+\use_package amssymb 1
+\use_package cancel 1
+\use_package esint 1
+\use_package mathdots 1
+\use_package mathtools 1
+\use_package mhchem 1
+\use_package stackrel 1
+\use_package stmaryrd 1
+\use_package undertilde 1
+\cite_engine basic
+\cite_engine_type default
+\biblio_style plain
+\use_bibtopic false
+\use_indices false
+\paperorientation portrait
+\suppress_date false
+\justification true
+\use_refstyle 1
+\index Index
+\shortcut idx
+\color #008000
+\end_index
+\secnumdepth 3
+\tocdepth 3
+\paragraph_separation indent
+\paragraph_indentation default
+\quotes_language english
+\papercolumns 1
+\papersides 1
+\paperpagestyle default
+\tracking_changes false
+\output_changes false
+\html_math_output 0
+\html_css_as_file 0
+\html_be_strict false
+\end_header
+
+\begin_body
+
+\begin_layout Title
+Boundary layer flow
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\begin{equation}
+\frac{\partial u}{\partial t}-\frac{1}{s^{2}}\frac{\partial^{2}u}{\partial y^{2}}=K(t)
+\end{equation}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula $y=0:u=0$
+\end_inset
+
+, 
+\begin_inset Formula $y=1,\frac{\partial u}{\partial y}=0$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Discretization, FTCD:
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\begin{equation}
+\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}-\frac{1}{s^{2}}\frac{u_{i+1}^{n}-2u_{i}^{n}-u_{i-1}^{n}}{\Delta y^{2}}=K^{n}
+\end{equation}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Rewriting:
+\begin_inset Note Note
+status open
+
+\begin_layout Plain Layout
+\begin_inset Formula $u_{i}^{n+1}-u_{i}^{n}-\frac{\Delta t}{\Delta y^{2}s^{2}}\frac{u_{i+1}^{n}-2u_{i}^{n}-u_{i-1}^{n}}{}=K^{n}\Delta t$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $u_{i}^{n+1}=\Delta tK^{n}+u_{i}^{n}+\frac{\Delta t}{\Delta y^{2}s^{2}}\left(u_{i+1}^{n}-2u_{i}^{n}-u_{i-1}^{n}\right)$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $u_{0}=0$
+\end_inset
+
+
+\end_layout
+
+\begin_layout Plain Layout
+and
+\end_layout
+
+\begin_layout Plain Layout
+\begin_inset Formula $u_{N}-u_{n-1}=0$
+\end_inset
+
+
+\end_layout
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+\begin_inset Formula 
+\[
+\]
+
+\end_inset
+
+
+\end_layout
+
+\end_body
+\end_document

blflow.py

@@ -0,0 +1,73 @@
+#!/usr/bin/python
+
+# Boundary layer flow
+from numpy import *
+import time 
+import matplotlib
+matplotlib.use('TkAgg')
+# from matplotlib.pylab import *
+import pylab as p
+
+# import matplotlib.animation as animation
+def K(t):                                 #Forcing function
+    return (1-exp(-0.1*t))*cos(t)
+
+
+s=10
+#Define domain
+n=50                                      #Number of gridpoints
+y=linspace(0,1,n)
+
+dy=y[1]-y[0]
+dt=0.0005
+
+l=(dt/(s**2*dy**2))
+
+hnu=exp(-sqrt(1j)*s*y)
+# fnu=(1-1j)/s
+fnu=0
+
+def u_ex(tn):
+    return (((1-hnu)/(1-fnu))*exp(1j*(tn))/1j).real
+    
+
+def u_np1(un,tn,dt):
+    Kn=K(tn)
+    unp1=un
+    unp1[0]=0                             #Velocity zero ver here
+    for i in range(1,un.size-1):
+        unp1[i]=dt*Kn+un[i]+l*(un[i-1]-2*un[i]+un[i+1])
+    unp1[-1]=unp1[-2]                     #Approximate 'infinity' bc
+    return unp1
+
+un0=zeros(n,float)
+t=0
+un=un0
+# un.append(un0)
+
+# Make the plot
+p.ion()
+linefd, = p.plot(un0,y)
+linee, = p.plot(un0,y)
+p.legend(('Finite difference','Periodic exact'))
+p.ylim(0,1)
+p.xlim(-1.5,1.5)
+p.ylabel('y')
+p.xlabel('u')
+p.grid('on')
+i=0
+uold=un
+while(True):
+    t+=dt
+    uold=un
+    un=u_np1(uold,t,dt)
+    if(i%20==0):
+        linefd.set_xdata(un)
+        linee.set_xdata(u_ex(t))
+        p.draw()
+        # print("Time:",t)
+    i+=1
+
+
+
+