Ported "Spectral methods in MATLAB" to python. Ported Series solution…

…s of Laplace's equation to python
mutablehurdle · Jan 28, 2022 · 57a9cdc · 57a9cdc
1 parent e7cf1cd
commit 57a9cdc
Show file tree

Hide file tree

Showing 60 changed files with 3,193 additions and 0 deletions.
diff --git a/series-solution-laplace/.ipynb_checkpoints/LaplaceManyDisks-checkpoint.ipynb b/series-solution-laplace/.ipynb_checkpoints/LaplaceManyDisks-checkpoint.ipynb
diff --git a/series-solution-laplace/LaplaceDisk.ipynb b/series-solution-laplace/LaplaceDisk.ipynb
diff --git a/series-solution-laplace/LaplaceDisk.py b/series-solution-laplace/LaplaceDisk.py
@@ -0,0 +1,53 @@
+
+
+
+#import necessary modules
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+#define center of disk using complex coordinates (real numbers are x, imaginary are y)
+
+c = 3 +1j #center of disk
+r = 1    #radius of disk
+
+#use power series expansion to solve Laplace's equation
+
+N =10      #number of expansion terms
+npts = 3*N #number of sample points
+
+npts_range=np.linspace(0,npts,npts)
+z = c + r*np.exp(2j*np.pi*npts_range/npts) #sample points along disk boundary
+
+rhs = -np.log(np.abs(z))+np.log(np.abs(z-c)) #right-hand side
+A = np.ones((npts,2*N+1))
+
+for k in np.arange(1,N+1):
+    A[:,2*k-1] = np.real((z-c)**(-k))
+    A[:,2*k] = np.imag((z-c)**(-k))
+
+a = np.linalg.lstsq(A,rhs)[0]
+
+
+#plot disk to check that it is centered where you want it, with the radius given and sample points
+plt.plot(np.real(z-c),np.imag(z-c),'.')
+
+
+#Define power series solution to Laplace's equation in complex numbers
+def disk1fun(z):
+    u = np.log(np.abs(z)) - np.log(np.abs(z-c)) + a[0]
+
+    for k in np.arange(1,N+1):
+        u+=a[2*k-1]*np.real((z-c)**(-k))+a[2*k]*np.imag((z-c)**(-k))
+    u[abs(z-c)<=r]=np.NaN
+
+    return u
+
+
+#contour plot the solution on a grid
+xx,yy = np.linspace(-5,5,145),np.linspace(-4,4,115)
+[xx,yy] = np.meshgrid(xx,yy)
+zz = xx + 1j*yy
+uu = disk1fun(zz)
+levels = np.linspace(-3,-.25,10)
+plt.contourf(xx,yy,uu)
diff --git a/series-solution-laplace/LaplaceManyDisks.ipynb b/series-solution-laplace/LaplaceManyDisks.ipynb
diff --git a/...ion-laplace/SOLVING LAPLACE PROBLEMS WITH CORNER SINGULARITIES VIA RATIONAL FUNCTIONS.pdf b/...ion-laplace/SOLVING LAPLACE PROBLEMS WITH CORNER SINGULARITIES VIA RATIONAL FUNCTIONS.pdf
diff --git a/series-solution-laplace/series solutions of laplace problems.pdf b/series-solution-laplace/series solutions of laplace problems.pdf
diff --git a/spectral-methods-python/.ipynb_checkpoints/p16_notebook-checkpoint.ipynb b/spectral-methods-python/.ipynb_checkpoints/p16_notebook-checkpoint.ipynb
@@ -0,0 +1,123 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAZ8ElEQVR4nO2dXawdV3mGn6/5pUBxEg6WZTsNCEsoFyXER2kQqKKJQMRFJBcBBaHGQpYstakEohJNWqkVUi+gFwQiVVCrQTUVP0n5UawoLbhOUNULAjb5x6Q5RIliK4lNSAIVom3g68VeE28O53jW7LP2nnfNfI+0dWavPX73mpm13nnXmtljc3eCIBgvv9V3BYIg6JcwgSAYOWECQTBywgSCYOSECQTByAkTCIKR07sJmNm7zexRM1sxsxt7rMfnzeyEmT08VXa+mR00s8fS3/NSuZnZLanOD5rZpQuq43Yzu8fMfmBmj5jZh9XqaWbnmtl3zeyBVMePp/LXm9m9qS63mdnZqfyc9H4lfX7RvOs4VdczzOw+M7tTsY5m9oSZPWRm95vZ4VRW/li7e28v4AzgR8AbgLOBB4CLe6rLHwCXAg9Plf0dcGNavhH4ZFreBfwrYMDlwL0LquMW4NK0/Grgv4CLleqZvutVafks4N703bcD16XyzwF/kpb/FPhcWr4OuG2Bx/yjwJeAO9N7qToCTwCvXVVW/FgvZGefZiPfCnxz6v1NwE091ueiVSbwKLAlLW8BHk3L/wB8YK31FlzfO4B3qtYT+G3g+8DvAz8Gzlx93IFvAm9Ny2em9WwBddsGHAKuAO5MnUetjmuZQPFj3fdwYCvw1NT7Y6lMhc3u/nRafgbYnJZ7r3eKpG9hcqaVqmeK2fcDJ4CDTNLeC+7+0hr1eLmO6fMXgQvmXUfg08DHgF+l9xcI1tGBb5nZETPbm8qKH+szS9R0DLi7m5nEPdZm9irga8BH3P2nZvbyZwr1dPdfApeY2SbgG8Cb+qzPaszsPcAJdz9iZu/ouz6n4e3uftzMXgccNLMfTn9Y6lj3nQSOA9un3m9LZSo8a2ZbANLfE6m8t3qb2VlMDOCL7v511XoCuPsLwD1MovUmM2tOOtP1eLmO6fPXAM/NuWpvA95rZk8AX2EyJPiMWB1x9+Pp7wkmZnoZczjWfZvA94AdaVb2bCaTLgd6rtM0B4DdaXk3kzF4U359mpG9HHhxKqLNDZuc8m8Fjrr7pxTraWZLKQFgZq9gMmdxlIkZXLtOHZu6Xwvc7WlQOy/c/SZ33+buFzFpc3e7+weV6mhmrzSzVzfLwLuAh5nHsV7EBFHL5McuJrPcPwL+qsd6fBl4Gvg/JuOpPUzGfYeAx4B/B85P6xrw96nODwHLC6rj25mMEx8E7k+vXUr1BH4PuC/V8WHgr1P5G4DvAivAvwDnpPJz0/uV9PkbFnzc38GpqwMydUx1eSC9Hmn6xjyOtSWBIAhGSt/DgSAIeiZMIAhGTphAEIycMIEgGDlzMQET+VFQEATtFDcBMzuDyaWKq5j8uOUDZnZxy7/Ze7rPFYg6lqOGeo6pjvNIApcBK+7+uLv/L5M7sq5u+TfyO5yoY0lqqOdo6jgPE+j9xzVBEOTT2w+IUpRpnGxn2w8hdu7cma195MiRTuvnaF144YUsLy9v6M6qkvVaixJ1XI/SdS9Vz3kca9j4vpxXvabpWscjR4782N2XVpfPwwSyfsjg7vuAfQA5v4Q6fPhwp0qYGaXuhgyt0BoCZvbkWuXzGA5I/CjI3Zn+iW1otWuVog+tnP1QUisXVa1pipuATx668GdMnsZyFLjd3R8p/T2ZdalaKxp2Psrmqqg1jcQPiHKGAxupp2qkDC1tVPfDrFpmdsTdl1eXj+KOQVU3HkO6qDmpKLebkvtrFCYA2gdhkVpdjGcMWm0oa5VqN6MxAdA1gtrTRUmtoaeLklql2o2ECezcubP96SdmWa82lDuvYkNbtFaki8VrSZjAolHtvMpapYh0oac1ShOAaNhdtXJR1FI2VwWt0ZpAg2KjzdEaQmyuOYIrdN5SWqM3AYWDsJ5WKVTTRc1JRbnddNUavQlA/wfhdES6GE+6KKnVpd2ECSRUjaD2dFFSa+jpoqRWl3YTJjCFcudVbGiL1hpCulBMKmECq1DtvGPRKoVqulBMKmECaxCxuT+tXFRNrEYtmV8RlvyFVRtdGm3fv/wKrfpQ3Q+m/CvCnTt3Fj1LlELZ2UuhmggiXZTXWg8JE4DFbOwsKB/QmJSrf1Ju0VprIWMCEEbQt1YpxpBUclHVmkbKBGD+rjcryp1XtaFFUtFNKtPImUBDJILQWq1VCtV00UdSAWETGMvQoBSqjbHmCK5siCX3l6wJwDiMACI296nVhrJWqXYjbQIwDiOoPV2U1Bp6uiipVard9PbfkHWh2dicg5+zTu6Oa9PqUq+c7xK+yWRhWl0TgaJWG2pa8kmgIRJBv1qliHShp1WNCUBcPuxTC3QabVetIcxdzHMepCoTaIhEEFqrtUqhmi7mmVSqNIGxDA1KodoYa47gyobYVatKE4BxGAFEbO5Tqw1lrS7tploTgHEYQe3poqTW0NNFSa0u7aZqE4AwglmoVWsI6UIxqVRvAhBGMCStUqimC8WkMggTAG0jWLRWNOx8lA1xUVqDfLxYSa02unQA1W0cutaiUd0PNuvjxczs82Z2wswenio738wOmtlj6e95qdzM7BYzWzGzB83s0pzKlXy8mHIiiNjcj1Yuim1wEe05ZzjwT8C7V5XdCBxy9x3AofQe4CpgR3rtBT6bW5Hadtws1DzBN4SJNMVJuUVrrUWrCbj7fwA/WVV8NbA/Le8Hrpkq/4JP+A6wycy25FYmjKBfrVKMIankoqo1zawTg5vd/em0/AywOS1vBZ6aWu9YKvsNzGyvmR02s8MnT558ubzPmyYWhaoRgG6jVU0Eqlpd2PDVgTSj17l27r7P3ZfdfXlpaek3PlftJKVQNQJlrVKopos+kgrMbgLPNjE//T2Ryo8D26fW25bKOqPaGEsSDbubVi6K7UZVC2Y3gQPA7rS8G7hjqvz6dJXgcuDFqWFDZ5R3XEkiNuvGZmWtUu2m9clCZvZl4B3Aa83sGPA3wCeA281sD/Ak8P60+l3ALmAF+DnwoY1WsNnYEjuwpFZJSm9jKfpKKm3rxr0ZE0q1m1YTcPcPrPPRlWus68ANG6rR2nXoZcflrJPrxjkNW7GhLVqrayJQ1GpDTaua24bHMDRQ3Ua1RttVa+hzFxvVqsYEQHdMVRJVIwCdRttVawhzF/OcB6nKBBpUO0kpVI1AWasUqulinkmlShNQbYwliYbdTSsXxXbTt1aVJgD977hFEbFZJzbXpNWl3VRrAjAOI6g9NpfUGnq6KKnVpd1UbQIQRjALtWoNIV0oJpXqTQDCCIakVQrVdKGYVAZhAhCXD4eiBZEuFq0VjxdboFYbcTvsfLQWjep+sFkfL7YIVB8vppwIFq0VsTmf2tqghAmA7o5TNQKI2NynVhuqWmshYwKg23lVjaD2SbmSWkNPF6W1ppEyAdCd4BuDEYBuo1VNBKpaXZAzgQbFzjsGI1DWKoVquugjqYCwCdQeTxeNcueNdKGrBcIm0FBrA+oD1c6rrFWKmpOKvAkoNyBVI6hZSzU2K2qVaoOtjxdToNnYEo2yL62cdXIPaK6W6A0rp9XqGnUVtdpQ05JPAg2qZ3HlRBCxeUKki9NTjQlAXD7sSs2TciW1VCflVC4fVmUCDYqddwxGoKxVCtV0Mc+kImkCcXtnWZQ7b6SL/hOBpAmMeXw2L1Q7r7JWKVSTSoOkCeSi3IBUjaBmrb5jc01aXdpg1SYAup1X1QhAs9HmaA0hgisOdas3AdDtvKpGoLqNY0gXikllECYAcfmwK8rbGOlisVrxeLFKtdrocsZR3UaFtjkLqvvB1B8vlovimGrRWiVRTQSRLsprrYeECTTUOqbqS6sUqp1XWasUClpSJqB80BW1SqLQGDeiFSeQ2Wk1ATPbbmb3mNkPzOwRM/twKj/fzA6a2WPp73mp3MzsFjNbMbMHzezSLhUawwSfqhGAbqOtNYKrXj6cJicJvAT8ubtfDFwO3GBmFwM3AofcfQdwKL0HuArYkV57gc/OUjHFDqeqVRLVbRxDuugjqUCGCbj70+7+/bT8M+AosBW4GtifVtsPXJOWrwa+4BO+A2wysy1dK1Z7PM3VCiPoRqSL8pOFneYEzOwi4C3AvcBmd386ffQMsDktbwWemvpnx1LZaq29ZnbYzA6fPHnydN/ZpYqnRVErjKBfrVLUnFSyTcDMXgV8DfiIu/901Rc50GkvuPs+d1929+WlpaXTrZetGZcPy6LceRUNfdFapfpGlgmY2VlMDOCL7v71VPxsE/PT3xOp/Diwfeqfb0tlG0J1TNWHlru3vsws65XzXYqdV1mrFIvSyrk6YMCtwFF3/9TURweA3Wl5N3DHVPn16SrB5cCLU8OGmVE+6IpaJamxYU8TJ5DTk/Og0bcBfww8ZGb3p7K/BD4B3G5me4Angfenz+4CdgErwM+BD81cu1XMY1KuhKaqVmlK1muRWl1j89C1VtNqAu7+n8B6NnPlGus7cMPMNcpAscOpapVEdRv7Shdt66pqrUbqjsGGmODTHhqoDn9qnJTrU6tB0gTGPD6bJoygX61SqM6DNEiaQIOqg9aaVEqi3HkVDX3RWl3ajbQJqM7+1p5USqHaeceiVQppE8hF+UApapVEtWGrmr7iCWQQJgC6k3KqWiWpOYIPYUi5Ua14vFhota4zz8tTi9BaNKr7weLxYt1R1hp6Ioh0UV5rPSRMoKHWMVUfWmEE/WqVQkFLygRUHTQuH3ZDufMqGnrfWlIm0AXVGdvak0opVDvvWLS6IGsCqjt36FolqT2Cj+EEAsImoBqbVSf4xmAEEBN8pbVA2AQaam1AuYQR9KtVipqTirwJjGFSTjVdlETVCEDz5JCjVapv5DxURIKcmyZqv6mllFbTSXK0ctbJbdg5x6ePbVy0VikWdXKQSQK1jqnGkFRKopoIxpIu1kLGBGoeU/WhlUsMDUKrDRkT6ILizh2DVklqn+Ab0glE0gRqjc2qE3xjMAKoc1KuT60GSROoPTYrdt4xGEHt6aKkVpf9IGkCDaoOWnNSCSPoRq1aXdqNtAmojqlqTiphBMPRKoW0CeSifKCGrlUS1disqpVLm9YgTAB0Y/MYtEqjGpuHqhWPFwutolpt9HFX56IRPj7xeLGuKGsNPRHUPMGnmgjWQ8IEGmodU/WhFUbQr1YpFLSkTEDVQePyYX+oGgFonhxm0ZIygS6oztjWnFTCCIaj1QVZE1DduaHVD3H5sJtWF2RNQDk2l9RS7LyqRgAxwVdaCzJMwMzONbPvmtkDZvaImX08lb/ezO41sxUzu83Mzk7l56T3K+nzi7JrPONG1Kyl2nlVjSAm+E5RKl3kJIH/Aa5w9zcDlwDvNrPLgU8CN7v7G4HngT1p/T3A86n85rTezIxlUq6kVhhBN2rVKtU3Wk3AJ/x3entWejlwBfDVVL4fuCYtX53ekz6/0grsGdUxlaJWX0bg7q0vM8t6lazXorVKsSitrDkBMzvDzO4HTgAHgR8BL7j7S2mVY8DWtLwVeCp98UvAi8AFa2juNbPDZnb45MmT1Y6pxpJUhp4IxpIu1iLLBNz9l+5+CbANuAx400a/2N33ufuyuy8vLS1VPabqQysX1U5SClUjUNVai05XB9z9BeAe4K3AJjNrnla8DTielo8D2wHS568BnitS24Tizg2t/ogTSDet1eRcHVgys01p+RXAO4GjTMzg2rTabuCOtHwgvSd9frd3PEo1x+aSWoqdV9UIoM5JuT61GnL+34EtwH4zO4OJadzu7nea2Q+Ar5jZ3wL3Abem9W8F/tnMVoCfANd1rVQX1yvV8RS1mg43ZK2SlN7GUvSVVHLXbTUBd38QeMsa5Y8zmR9YXf4L4H1Z395C24Z0dVBFrTbmkS7CCPKpVavL98jeMQi6Y6qY4CuvVRLVbVS9FCltAg2qYypVrTZUtUqi3HnVTg5VmEDM/nbTykVVqxSqnVdNKx4vFlqyWm0otN21EN6n8XixrihrKZ1J5qFVmrh8uD4SJtAQsTkf1c6ragSqk3IKWlImoOqgY5ngCyPoxlC0pEygC6qTcjUnlTCC4Wh1oVoTaFA9i6tqtaGqVRLlztvHyaF6E1BttJEu5qNVCtXO24eWvAmoNsaaD3rNWiWJE8gEeROI2NxdS7HzqhoBaJ4ccrRK9Y2cXxFKkHPTxDx+YVWjluoPhbpo5ayT2+Fy2s0Yfn243v6SSQKqE2mqWm2MIV2UZCyXD9dCxgRqHlP1oZWLYucdgxGoaq2FjAk01DimGopWG6paJVHtvPM0AjkTiNs7TzH0dFFaqxSqnXdeRiBnAhANu0G1AalqlWQMJ5AGSRNoUI3NNUdwxc6ragSgeXLI0erSbqRNQDU215xUVDuvqhGobmPJk4O0CTSonsVVtdpQ1goj6EYJrSpMQHVMpaqVi6JWGMHitaq5YxB0H9k8Bq02+tBatAkP4U7MtagiCTSojqnGkFRyUdUqhdpZvISWlAlEw85HpQHVolWS2k8gq5EygSFMytU8wafYeVWNADRPDrNoSZlAF1Rjc81JRbXzqhqB6jZ2PTlUawINqmdxVa02lLXCCLqRq1W9CSg32hyGni5KaoURzEerehNoUGy0OVpDSBe1JpWS1GwEgzGBuHx4ipgH6YdajaCam4Vybobo0hgXddNRH1pKN6IoaOWsk9vhSrXBEt/VVWu9bcxOAmZ2hpndZ2Z3pvevN7N7zWzFzG4zs7NT+Tnp/Ur6/KJM/daN6FDXwWu1oTopp6pVmppST5fhwIeBo1PvPwnc7O5vBJ4H9qTyPcDzqfzmtF4rEZu7aeWi2OFUtUpS0zZmmYCZbQP+CPjH9N6AK4CvplX2A9ek5avTe9LnV1qHLYhJufGki5JaYQSzk5sEPg18DPhVen8B8IK7v5TeHwO2puWtwFMA6fMX0/q/hpntNbPDZnb45MmTL5fHBN8php4uSmqFEcxOqwmY2XuAE+5+pOQXu/s+d1929+WlpaXVn0k2tEVrDSFd1JpUSqJuBDlXB94GvNfMdgHnAr8DfAbYZGZnprP9NuB4Wv84sB04ZmZnAq8BnutasSHMSudolaLklRFVrVxKapVCtQ1CRhJw95vcfZu7XwRcB9zt7h8E7gGuTavtBu5IywfSe9Lnd/uMtVWNzTVHcNWzkqpWSVSHuhu5WegvgI+a2QqTMf+tqfxW4IJU/lHgxo1VUTeeqmq1oawVRtCNElqdbhZy928D307LjwOXrbHOL4D3bbhmv665cK2Izf1oKcfmUqhtYzV3DEI8XqxPrTaUtUqbcM53lqpXDhvVquq3A6pjqrh8qK01hqHBRrSqMgHQHFP1oTWEuYta50FKomAE1ZkAaOy4RWiVQjVd1J5UStF3G6zSBEC3MdbcsPtujLVplaTPE0i1JtCgGk9VtdpQ1goj6EauVvUmoNxocxh6uiipFUYwH63qTaBBsdHmaA0hXdScVMIIBmQCtU/KldSKeZB8wgjAFC6dmNmsPy9YS0vuLrjQ6ldr0eR23pztK7xPj7j78upyiSSwc+dOSWdX1iqFaiKoOV2UZBFJRcIEIBp2V61cVE1s6FolmXe9ZEygISblxjMpV1IrjGB25Eyg9thcUmvo6aKkVhjB7MiZAMTvAxqGkC5qTipjMQJJEwBdZ689qcQ8SD5jMQJZEwDdxhgNO7T6puQJRNoEGlTjqapWG8paip1X1QigzMmhChOI2NxNKxdFLdXOq2oEJeoVjxdbh65n8aFrtaGsVard5Got+uSw0W2sIgk01D4pV1Ir0kU+kQhOT1UmAHH5sGEIcxc1z4MMyQiqMwHQdfbak0rMg+QzJCOo0gRAtzFGww6tvul6AqnWBBpU46mqVhvKWoqdV9UIIB4vVo3W0NNFSS3VzqtqBLn1qt4EGhQbbY7WENJFzUkljGBAJlD7pFxJrZgHySeMIB4vFloj0Fo0pfdDDjnfZ/F4sWFplUI1EUS6WBwSJgDRsLtq5aLYsENLCxkTaIhJufFMypXUUuy8tRhBlgmY2RNm9pCZ3W9mh1PZ+WZ20MweS3/PS+VmZreY2YqZPWhml3apUO2xuaTW0NNFSS3VzluDEXRJAn/o7pdMTSzcCBxy9x3AofQe4CpgR3rtBT7btVLx+4AJQ0gXNSeVsRjBRoYDVwP70/J+4Jqp8i/4hO8Am8xsS1dx1YNQe1KJeZB8VNtNaXJNwIFvmdkRM9ubyja7+9Np+Rlgc1reCjw19W+PpbLOqB6EmpPKEBJBzelCkdyHirzd3Y+b2euAg2b2w+kP3d3NrNMWJjPZC3DhhReuu17T4UrswLFolaKvpNK2bkmtXFS1SpCVBNz9ePp7AvgGcBnwbBPz098TafXjwPapf74tla3W3Ofuy+6+vLS01Pb9OdXMQjXq1hybS2qpJjZVrRK0moCZvdLMXt0sA+8CHgYOALvTaruBO9LyAeD6dJXgcuDFqWHDhojYHLF5Fi3Fzpur5e5ZLzNrfa1HznBgM/CNJHIm8CV3/zcz+x5wu5ntAZ4E3p/WvwvYBawAPwc+lPEdWURsPkXE5nyUh24KQ4NWE3D3x4E3r1H+HHDlGuUO3FCkdmvXR7KhLVqrayIYulYb80gXQzECuTsGc6g94uVqlUJ17qLmeRDVdjMLVZoA6B6EuHw4HK02VLW6Uq0JgG7nVdYqxRiSSi6qWrlUbQIQDburVi6KWsrmqqiVS/Um0KDYaHO0hhCba47gip130UYwGBOI2HyKmODLR7XzLtIIVB4v9jPg0b7r0cJrgR/3XYkWaqgj1FHPIdbxd939N27PVfkPSR/1NZ59poSZHY46lqGGeo6pjoMZDgRBMBthAkEwclRMYF/fFcgg6liOGuo5mjpKTAwGQdAfKkkgCIKeCBMIgpETJhAEIydMIAhGTphAEIyc/wegafmDpYEJkAAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "'''\n",
+    "Spectral methods in MATLAB.\n",
+    "compare to Trefethen p16.m\n",
+    "'''\n",
+    "\n",
+    "# Poisson equation on [-1,1]x[-1,1] with u=o on boundary\n",
+    "\n",
+    "from mayavi import mlab\n",
+    "from numpy import *\n",
+    "import time\n",
+    "from cheb import *\n",
+    "from numpy.linalg import matrix_power,solve\n",
+    "from matplotlib import pyplot as plt\n",
+    "from scipy.interpolate import interp2d,griddata,bisplev,bisplrep\n",
+    "\n",
+    "# Set up grids and tensor product Laplacian and solve for u:\n",
+    "\n",
+    "N = 24\n",
+    "D,x = cheb(N)\n",
+    "y = x\n",
+    "xx,yy = meshgrid(x[1:N],y[1:N])\n",
+    "xx = hstack(xx[:]); yy = hstack(yy[:]);                   # stretch 2D grids to 1D vectors\n",
+    "f = 10*sin(8*xx*(yy-1))\n",
+    "D2 = matrix_power(D, 2)\n",
+    "D2 = D2[1:N,1:N]\n",
+    "I = identity(N-1)\n",
+    "L = kron(I,D2) + kron(D2,I)          # Laplacian\n",
+    "plt.figure(1)\n",
+    "plt.spy(L)\n",
+    "tic = time.time(); u = solve(L,f); toc = time.time() - tic;      # Solve problem and watch the clock\n",
+    "\n",
+    "# Reshape long 1D results onto 2D grid:\n",
+    "uu = zeros((N+1,N+1))\n",
+    "uu[1:N,1:N]=u.reshape(N-1,N-1)\n",
+    "xx,yy = meshgrid(x,y)\n",
+    "value = uu[int(N/4.+1),int(N/4.+1)]\n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 47,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "uuu = interp2d(x, y, uu, kind='cubic')\n",
+    "\n",
+    "mlab.surf(uuu(arange(-1,1,0.04),arange(-1,1,0.04)),warp_scale=\"auto\")\n",
+    "mlab.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(50, 50)"
+      ]
+     },
+     "execution_count": 45,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "uuu(arange(-1,1,0.04),arange(-1,1,0.04)).shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/spectral-methods-python/.ipynb_checkpoints/p19_notebook-checkpoint.ipynb b/spectral-methods-python/.ipynb_checkpoints/p19_notebook-checkpoint.ipynb
@@ -0,0 +1,109 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<mayavi.modules.surface.Surface at 0x10c826e30>"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "'''\n",
+    "Spectral methods in MATLAB.\n",
+    "compare to Trefethen p18.m\n",
+    "'''\n",
+    "\n",
+    "# 2nd order wave eq. on Chebyshev grid\n",
+    "\n",
+    "from numpy import *\n",
+    "from matplotlib import pyplot as plt\n",
+    "from mpl_toolkits.mplot3d import axes3d\n",
+    "from chebfft import *\n",
+    "from mayavi import mlab\n",
+    "\n",
+    "N = 80\n",
+    "x = cos(pi*arange(0,N+1)/N)\n",
+    "dt = 8.0/N**2\n",
+    "v = exp(-200*x**2)\n",
+    "vold = exp(-200*(x-dt)**2)\n",
+    "tmax = 4\n",
+    "tplot = 0.075\n",
+    "plotgap = int(round(tplot/dt))\n",
+    "dt = tplot/plotgap\n",
+    "nplots = int(round(tmax/tplot))\n",
+    "plotdata = vstack((v, zeros((nplots,N+1))))\n",
+    "tdata = 0\n",
+    "for i in range(0,nplots):\n",
+    "    for n in range(0,plotgap):\n",
+    "        w = chebfft(chebfft(v)).T\n",
+    "        w[0] = 0\n",
+    "        w[N] = 0\n",
+    "        vnew = 2*v - vold +dt**2*w\n",
+    "        vold = v\n",
+    "        v = vnew\n",
+    "    plotdata[i+1,:] = v\n",
+    "    tdata = vstack((tdata, dt*i*plotgap))\n",
+    "\n",
+    "# Plot results\n",
+    "\n",
+    "# fig = plt.figure()\n",
+    "# ax = axes3d.Axes3D(fig)\n",
+    "X, Y = meshgrid(x, tdata)\n",
+    "\n",
+    "mlab.surf(plotdata,warp_scale=\"auto\")\n",
+    "\n",
+    "# ax.plot_wireframe(X,Y,plotdata)\n",
+    "# ax.set_xlim(-1, 1)\n",
+    "# ax.set_ylim(0, tmax)\n",
+    "# ax.set_zlim(-2, 2)\n",
+    "# plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mlab.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}