# Vibration modes of a membrane in parabolic coordinates



In [1]:

%matplotlib notebook




In [2]:

import numpy as np
from scipy.linalg import eigh
from sympy import (symbols, lambdify, init_printing,
expand, Matrix, diff, integrate)
from sympy.utilities.lambdify import lambdify
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D




In [3]:

x, y, r, s= symbols('x y r s')
init_printing()



## Boundary conditions

The boundary conditions are satisfied multiplying by $b(x, y)$.



In [4]:

b = lambda x, y: (2*y - x**2 + 1)*(2*y - 1 + x**2)
bound = b(x, y)




In [5]:

b_num = lambdify((x,y), bound, "numpy")
X, Y = np.mgrid[-1:1:200j, -0.5:0.5:200j]
bound_num = b_num(X,Y)




In [6]:

fig = plt.figure(figsize=(6, 3))
plt.contourf(X, Y, bound_num, 12, cmap="RdYlBu", vmin=-1, vmax=1)
plt.contour(X, Y, bound_num, [0], colors="black")
plt.axis("image");




var element = $('#7ef08960-9b93-4271-bee0-d0a306e8f464'); /* Put everything inside the global mpl namespace */ window.mpl = {}; mpl.get_websocket_type = function() { if (typeof(WebSocket) !== 'undefined') { return WebSocket; } else if (typeof(MozWebSocket) !== 'undefined') { return MozWebSocket; } else { alert('Your browser does not have WebSocket support. ' + 'Please try Chrome, Safari or Firefox ≥ 6. ' + 'Firefox 4 and 5 are also supported but you ' + 'have to enable WebSockets in about:config.'); }; } mpl.figure = function(figure_id, websocket, ondownload, parent_element) { this.id = figure_id; this.ws = websocket; this.supports_binary = (this.ws.binaryType != undefined); if (!this.supports_binary) { var warnings = document.getElementById("mpl-warnings"); if (warnings) { warnings.style.display = 'block'; warnings.textContent = ( "This browser does not support binary websocket messages. " + "Performance may be slow."); } } this.imageObj = new Image(); this.context = undefined; this.message = undefined; this.canvas = undefined; this.rubberband_canvas = undefined; this.rubberband_context = undefined; this.format_dropdown = undefined; this.image_mode = 'full'; this.root =$('<div/>');
this._root_extra_style(this.root)
this.root.attr('style', 'display: inline-block');

$(parent_element).append(this.root); this._init_header(this); this._init_canvas(this); this._init_toolbar(this); var fig = this; this.waiting = false; this.ws.onopen = function () { fig.send_message("supports_binary", {value: fig.supports_binary}); fig.send_message("send_image_mode", {}); if (mpl.ratio != 1) { fig.send_message("set_dpi_ratio", {'dpi_ratio': mpl.ratio}); } fig.send_message("refresh", {}); } this.imageObj.onload = function() { if (fig.image_mode == 'full') { // Full images could contain transparency (where diff images // almost always do), so we need to clear the canvas so that // there is no ghosting. fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height); } fig.context.drawImage(fig.imageObj, 0, 0); }; this.imageObj.onunload = function() { fig.ws.close(); } this.ws.onmessage = this._make_on_message_function(this); this.ondownload = ondownload; } mpl.figure.prototype._init_header = function() { var titlebar =$(
'<div class="ui-dialog-titlebar ui-widget-header ui-corner-all ' +
'ui-helper-clearfix"/>');
var titletext = ( '<div class="ui-dialog-title" style="width: 100%; ' + 'text-align: center; padding: 3px;"/>'); titlebar.append(titletext) this.root.append(titlebar); this.header = titletext[0]; } mpl.figure.prototype._canvas_extra_style = function(canvas_div) { } mpl.figure.prototype._root_extra_style = function(canvas_div) { } mpl.figure.prototype._init_canvas = function() { var fig = this; var canvas_div =('<div/>');

canvas_div.attr('style', 'position: relative; clear: both; outline: 0');

function canvas_keyboard_event(event) {
return fig.key_event(event, event['data']);
}

canvas_div.keydown('key_press', canvas_keyboard_event);
canvas_div.keyup('key_release', canvas_keyboard_event);
this.canvas_div = canvas_div
this._canvas_extra_style(canvas_div)
this.root.append(canvas_div);

var canvas = $('<canvas/>'); canvas.addClass('mpl-canvas'); canvas.attr('style', "left: 0; top: 0; z-index: 0; outline: 0") this.canvas = canvas[0]; this.context = canvas[0].getContext("2d"); var backingStore = this.context.backingStorePixelRatio || this.context.webkitBackingStorePixelRatio || this.context.mozBackingStorePixelRatio || this.context.msBackingStorePixelRatio || this.context.oBackingStorePixelRatio || this.context.backingStorePixelRatio || 1; mpl.ratio = (window.devicePixelRatio || 1) / backingStore; var rubberband =$('<canvas/>');
rubberband.attr('style', "position: absolute; left: 0; top: 0; z-index: 1;")

var pass_mouse_events = true;

canvas_div.resizable({
start: function(event, ui) {
pass_mouse_events = false;
},
resize: function(event, ui) {
fig.request_resize(ui.size.width, ui.size.height);
},
stop: function(event, ui) {
pass_mouse_events = true;
fig.request_resize(ui.size.width, ui.size.height);
},
});

function mouse_event_fn(event) {
if (pass_mouse_events)
return fig.mouse_event(event, event['data']);
}

rubberband.mousedown('button_press', mouse_event_fn);
rubberband.mouseup('button_release', mouse_event_fn);
// Throttle sequential mouse events to 1 every 20ms.
rubberband.mousemove('motion_notify', mouse_event_fn);

rubberband.mouseenter('figure_enter', mouse_event_fn);
rubberband.mouseleave('figure_leave', mouse_event_fn);

canvas_div.on("wheel", function (event) {
event = event.originalEvent;
event['data'] = 'scroll'
if (event.deltaY < 0) {
event.step = 1;
} else {
event.step = -1;
}
mouse_event_fn(event);
});

canvas_div.append(canvas);
canvas_div.append(rubberband);

this.rubberband = rubberband;
this.rubberband_canvas = rubberband[0];
this.rubberband_context = rubberband[0].getContext("2d");
this.rubberband_context.strokeStyle = "#000000";

this._resize_canvas = function(width, height) {
// Keep the size of the canvas, canvas container, and rubber band
// canvas in synch.
canvas_div.css('width', width)
canvas_div.css('height', height)

canvas.attr('width', width * mpl.ratio);
canvas.attr('height', height * mpl.ratio);
canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');

rubberband.attr('width', width);
rubberband.attr('height', height);
}

// Set the figure to an initial 600x600px, this will subsequently be updated
// upon first draw.
this._resize_canvas(600, 600);

// Disable right mouse context menu.
$(this.rubberband_canvas).bind("contextmenu",function(e){ return false; }); function set_focus () { canvas.focus(); canvas_div.focus(); } window.setTimeout(set_focus, 100); } mpl.figure.prototype._init_toolbar = function() { var fig = this; var nav_element =$('<div/>');
nav_element.attr('style', 'width: 100%');
this.root.append(nav_element);

// Define a callback function for later on.
function toolbar_event(event) {
return fig.toolbar_button_onclick(event['data']);
}
function toolbar_mouse_event(event) {
return fig.toolbar_button_onmouseover(event['data']);
}

for(var toolbar_ind in mpl.toolbar_items) {
var name = mpl.toolbar_items[toolbar_ind][0];
var tooltip = mpl.toolbar_items[toolbar_ind][1];
var image = mpl.toolbar_items[toolbar_ind][2];
var method_name = mpl.toolbar_items[toolbar_ind][3];

if (!name) {
// put a spacer in here.
continue;
}
var button = $('<button/>'); button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' + 'ui-button-icon-only'); button.attr('role', 'button'); button.attr('aria-disabled', 'false'); button.click(method_name, toolbar_event); button.mouseover(tooltip, toolbar_mouse_event); var icon_img =$('<span/>');

var tooltip_span = $('<span/>'); tooltip_span.addClass('ui-button-text'); tooltip_span.html(tooltip); button.append(icon_img); button.append(tooltip_span); nav_element.append(button); } var fmt_picker_span =$('<span/>');

var fmt_picker = $('<select/>'); fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content'); fmt_picker_span.append(fmt_picker); nav_element.append(fmt_picker_span); this.format_dropdown = fmt_picker[0]; for (var ind in mpl.extensions) { var fmt = mpl.extensions[ind]; var option =$(
'<option/>', {selected: fmt === mpl.default_extension}).html(fmt);
fmt_picker.append(option);
}

// Add hover states to the ui-buttons
$( ".ui-button" ).hover( function() {$(this).addClass("ui-state-hover");},
function() { $(this).removeClass("ui-state-hover");} ); var status_bar =$('<span class="mpl-message"/>');
nav_element.append(status_bar);
this.message = status_bar[0];
}

mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {
// Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,
// which will in turn request a refresh of the image.
this.send_message('resize', {'width': x_pixels, 'height': y_pixels});
}

mpl.figure.prototype.send_message = function(type, properties) {
properties['type'] = type;
properties['figure_id'] = this.id;
this.ws.send(JSON.stringify(properties));
}

mpl.figure.prototype.send_draw_message = function() {
if (!this.waiting) {
this.waiting = true;
this.ws.send(JSON.stringify({type: "draw", figure_id: this.id}));
}
}

mpl.figure.prototype.handle_save = function(fig, msg) {
var format_dropdown = fig.format_dropdown;
var format = format_dropdown.options[format_dropdown.selectedIndex].value;
}

mpl.figure.prototype.handle_resize = function(fig, msg) {
var size = msg['size'];
if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {
fig._resize_canvas(size[0], size[1]);
fig.send_message("refresh", {});
};
}

mpl.figure.prototype.handle_rubberband = function(fig, msg) {
var x0 = msg['x0'] / mpl.ratio;
var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;
var x1 = msg['x1'] / mpl.ratio;
var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;
x0 = Math.floor(x0) + 0.5;
y0 = Math.floor(y0) + 0.5;
x1 = Math.floor(x1) + 0.5;
y1 = Math.floor(y1) + 0.5;
var min_x = Math.min(x0, x1);
var min_y = Math.min(y0, y1);
var width = Math.abs(x1 - x0);
var height = Math.abs(y1 - y0);

fig.rubberband_context.clearRect(
0, 0, fig.canvas.width / mpl.ratio, fig.canvas.height / mpl.ratio);

fig.rubberband_context.strokeRect(min_x, min_y, width, height);
}

mpl.figure.prototype.handle_figure_label = function(fig, msg) {
}

mpl.figure.prototype.handle_cursor = function(fig, msg) {
var cursor = msg['cursor'];
switch(cursor)
{
case 0:
cursor = 'pointer';
break;
case 1:
cursor = 'default';
break;
case 2:
cursor = 'crosshair';
break;
case 3:
cursor = 'move';
break;
}
fig.rubberband_canvas.style.cursor = cursor;
}

mpl.figure.prototype.handle_message = function(fig, msg) {
fig.message.textContent = msg['message'];
}

mpl.figure.prototype.handle_draw = function(fig, msg) {
// Request the server to send over a new figure.
fig.send_draw_message();
}

mpl.figure.prototype.handle_image_mode = function(fig, msg) {
fig.image_mode = msg['mode'];
}

mpl.figure.prototype.updated_canvas_event = function() {
// Called whenever the canvas gets updated.
this.send_message("ack", {});
}

// A function to construct a web socket function for onmessage handling.
// Called in the figure constructor.
mpl.figure.prototype._make_on_message_function = function(fig) {
return function socket_on_message(evt) {
if (evt.data instanceof Blob) {
/* FIXME: We get "Resource interpreted as Image but
* transferred with MIME type text/plain:" errors on
* Chrome.  But how to set the MIME type?  It doesn't seem
* to be part of the websocket stream */
evt.data.type = "image/png";

/* Free the memory for the previous frames */
if (fig.imageObj.src) {
(window.URL || window.webkitURL).revokeObjectURL(
fig.imageObj.src);
}

fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(
evt.data);
fig.updated_canvas_event();
fig.waiting = false;
return;
}
else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == "data:image/png;base64") {
fig.imageObj.src = evt.data;
fig.updated_canvas_event();
fig.waiting = false;
return;
}

var msg = JSON.parse(evt.data);
var msg_type = msg['type'];

// Call the  "handle_{type}" callback, which takes
// the figure and JSON message as its only arguments.
try {
var callback = fig["handle_" + msg_type];
} catch (e) {
console.log("No handler for the '" + msg_type + "' message type: ", msg);
return;
}

if (callback) {
try {
// console.log("Handling '" + msg_type + "' message: ", msg);
callback(fig, msg);
} catch (e) {
console.log("Exception inside the 'handler_" + msg_type + "' callback:", e, e.stack, msg);
}
}
};
}

// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas
mpl.findpos = function(e) {
//this section is from http://www.quirksmode.org/js/events_properties.html
var targ;
if (!e)
e = window.event;
if (e.target)
targ = e.target;
else if (e.srcElement)
targ = e.srcElement;
if (targ.nodeType == 3) // defeat Safari bug
targ = targ.parentNode;

// jQuery normalizes the pageX and pageY
// pageX,Y are the mouse positions relative to the document
// offset() returns the position of the element relative to the document
var x = e.pageX - $(targ).offset().left; var y = e.pageY -$(targ).offset().top;

return {"x": x, "y": y};
};

/*
* return a copy of an object with only non-object keys
* we need this to avoid circular references
* http://stackoverflow.com/a/24161582/3208463
*/
function simpleKeys (original) {
return Object.keys(original).reduce(function (obj, key) {
if (typeof original[key] !== 'object')
obj[key] = original[key]
return obj;
}, {});
}

mpl.figure.prototype.mouse_event = function(event, name) {
var canvas_pos = mpl.findpos(event)

if (name === 'button_press')
{
this.canvas.focus();
this.canvas_div.focus();
}

var x = canvas_pos.x * mpl.ratio;
var y = canvas_pos.y * mpl.ratio;

this.send_message(name, {x: x, y: y, button: event.button,
step: event.step,
guiEvent: simpleKeys(event)});

/* This prevents the web browser from automatically changing to
* the text insertion cursor when the button is pressed.  We want
* to control all of the cursor setting manually through the
* 'cursor' event from matplotlib */
event.preventDefault();
return false;
}

mpl.figure.prototype._key_event_extra = function(event, name) {
// Handle any extra behaviour associated with a key event
}

mpl.figure.prototype.key_event = function(event, name) {

// Prevent repeat events
if (name == 'key_press')
{
if (event.which === this._key)
return;
else
this._key = event.which;
}
if (name == 'key_release')
this._key = null;

var value = '';
if (event.ctrlKey && event.which != 17)
value += "ctrl+";
if (event.altKey && event.which != 18)
value += "alt+";
if (event.shiftKey && event.which != 16)
value += "shift+";

value += 'k';
value += event.which.toString();

this._key_event_extra(event, name);

this.send_message(name, {key: value,
guiEvent: simpleKeys(event)});
return false;
}

mpl.figure.prototype.toolbar_button_onclick = function(name) {
this.handle_save(this, null);
} else {
this.send_message("toolbar_button", {name: name});
}
};

mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {
this.message.textContent = tooltip;
};
mpl.toolbar_items = [["Home", "Reset original view", "fa fa-home icon-home", "home"], ["Back", "Back to previous view", "fa fa-arrow-left icon-arrow-left", "back"], ["Forward", "Forward to next view", "fa fa-arrow-right icon-arrow-right", "forward"], ["", "", "", ""], ["Pan", "Pan axes with left mouse, zoom with right", "fa fa-arrows icon-move", "pan"], ["Zoom", "Zoom to rectangle", "fa fa-square-o icon-check-empty", "zoom"], ["", "", "", ""], ["Download", "Download plot", "fa fa-floppy-o icon-save", "download"]];

mpl.extensions = ["eps", "jpeg", "pdf", "png", "ps", "raw", "svg", "tif"];

mpl.default_extension = "png";var comm_websocket_adapter = function(comm) {
// Create a "websocket"-like object which calls the given IPython comm
// object with the appropriate methods. Currently this is a non binary
// socket, so there is still some room for performance tuning.
var ws = {};

ws.close = function() {
comm.close()
};
ws.send = function(m) {
//console.log('sending', m);
comm.send(m);
};
// Register the callback with on_msg.
comm.on_msg(function(msg) {
//console.log('receiving', msg['content']['data'], msg);
// Pass the mpl event to the overridden (by mpl) onmessage function.
ws.onmessage(msg['content']['data'])
});
return ws;
}

mpl.mpl_figure_comm = function(comm, msg) {
// This is the function which gets called when the mpl process
// starts-up an IPython Comm through the "matplotlib" channel.

var id = msg.content.data.id;
// Get hold of the div created by the display call when the Comm
// socket was opened in Python.
var element = $("#" + id); var ws_proxy = comm_websocket_adapter(comm) function ondownload(figure, format) { window.open(figure.imageObj.src); } var fig = new mpl.figure(id, ws_proxy, ondownload, element.get(0)); // Call onopen now - mpl needs it, as it is assuming we've passed it a real // web socket which is closed, not our websocket->open comm proxy. ws_proxy.onopen(); fig.parent_element = element.get(0); fig.cell_info = mpl.find_output_cell("<div id='" + id + "'></div>"); if (!fig.cell_info) { console.error("Failed to find cell for figure", id, fig); return; } var output_index = fig.cell_info[2] var cell = fig.cell_info[0]; }; mpl.figure.prototype.handle_close = function(fig, msg) { var width = fig.canvas.width/mpl.ratio fig.root.unbind('remove') // Update the output cell to use the data from the current canvas. fig.push_to_output(); var dataURL = fig.canvas.toDataURL(); // Re-enable the keyboard manager in IPython - without this line, in FF, // the notebook keyboard shortcuts fail. IPython.keyboard_manager.enable()$(fig.parent_element).html('<amp-img layout="responsive" width="500" height="300" src="' + dataURL + '" width="' + width + '">');
fig.close_ws(fig, msg);
}

mpl.figure.prototype.close_ws = function(fig, msg){
fig.send_message('closing', msg);
// fig.ws.close()
}

mpl.figure.prototype.push_to_output = function(remove_interactive) {
// Turn the data on the canvas into data in the output cell.
var width = this.canvas.width/mpl.ratio
var dataURL = this.canvas.toDataURL();
this.cell_info[1]['text/html'] = '<amp-img layout="responsive" width="500" height="300" src="' + dataURL + '" width="' + width + '">';
}

mpl.figure.prototype.updated_canvas_event = function() {
// Tell IPython that the notebook contents must change.
IPython.notebook.set_dirty(true);
this.send_message("ack", {});
var fig = this;
// Wait a second, then push the new image to the DOM so
// that it is saved nicely (might be nice to debounce this).
setTimeout(function () { fig.push_to_output() }, 1000);
}

mpl.figure.prototype._init_toolbar = function() {
var fig = this;

var nav_element = $('<div/>'); nav_element.attr('style', 'width: 100%'); this.root.append(nav_element); // Define a callback function for later on. function toolbar_event(event) { return fig.toolbar_button_onclick(event['data']); } function toolbar_mouse_event(event) { return fig.toolbar_button_onmouseover(event['data']); } for(var toolbar_ind in mpl.toolbar_items){ var name = mpl.toolbar_items[toolbar_ind][0]; var tooltip = mpl.toolbar_items[toolbar_ind][1]; var image = mpl.toolbar_items[toolbar_ind][2]; var method_name = mpl.toolbar_items[toolbar_ind][3]; if (!name) { continue; }; var button =$('<button class="btn btn-default" href="#" title="' + name + '"><i class="fa ' + image + ' fa-lg"></i></button>');
button.click(method_name, toolbar_event);
button.mouseover(tooltip, toolbar_mouse_event);
nav_element.append(button);
}

var status_bar = ('<span class="mpl-message" style="text-align:right; float: right;"/>'); nav_element.append(status_bar); this.message = status_bar[0]; // Add the close button to the window. var buttongrp =('<div class="btn-group inline pull-right"></div>');
var button = $('<button class="btn btn-mini btn-primary" href="#" title="Stop Interaction"><i class="fa fa-power-off icon-remove icon-large"></i></button>'); button.click(function (evt) { fig.handle_close(fig, {}); } ); button.mouseover('Stop Interaction', toolbar_mouse_event); buttongrp.append(button); var titlebar = this.root.find($('.ui-dialog-titlebar'));
titlebar.prepend(buttongrp);
}

mpl.figure.prototype._root_extra_style = function(el){
var fig = this
el.on("remove", function(){
fig.close_ws(fig, {});
});
}

mpl.figure.prototype._canvas_extra_style = function(el){
// this is important to make the div 'focusable
el.attr('tabindex', 0)
// reach out to IPython and tell the keyboard manager to turn it's self
// off when our div gets focus

// location in version 3
if (IPython.notebook.keyboard_manager) {
IPython.notebook.keyboard_manager.register_events(el);
}
else {
// location in version 2
IPython.keyboard_manager.register_events(el);
}

}

mpl.figure.prototype._key_event_extra = function(event, name) {
var manager = IPython.notebook.keyboard_manager;
if (!manager)
manager = IPython.keyboard_manager;

// Check for shift+enter
if (event.shiftKey && event.which == 13) {
this.canvas_div.blur();
// select the cell after this one
var index = IPython.notebook.find_cell_index(this.cell_info[0]);
IPython.notebook.select(index + 1);
}
}

mpl.figure.prototype.handle_save = function(fig, msg) {
}

mpl.find_output_cell = function(html_output) {
// Return the cell and output element which can be found *uniquely* in the notebook.
// Note - this is a bit hacky, but it is done because the "notebook_saving.Notebook"
// IPython event is triggered only after the cells have been serialised, which for
// our purposes (turning an active figure into a static one), is too late.
var cells = IPython.notebook.get_cells();
var ncells = cells.length;
for (var i=0; i<ncells; i++) {
var cell = cells[i];
if (cell.cell_type === 'code'){
for (var j=0; j<cell.output_area.outputs.length; j++) {
var data = cell.output_area.outputs[j];
if (data.data) {
// IPython >= 3 moved mimebundle to data attribute of output
data = data.data;
}
if (data['text/html'] == html_output) {
return [cell, data, j];
}
}
}
}
}

// Register the function which deals with the matplotlib target/channel.
// The kernel may be null if the page has been refreshed.
if (IPython.notebook.kernel != null) {
IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);
}



## Approximating functions



In [7]:

def w_fun(x, y, m, n):
""" Trial function. """
c = symbols('c:%d' % (m*n)) # This is the way of define the coefficients c_i
w = []
for i in range(0, m):
for j in range(0, n):
w.append(x**i * y**j)

return w, c

def u_fun(x, y, m, n):
""" Complete function. Contains the boundary and trial functions. """
w, c = w_fun(x, y, m, n)
return [b(x, y) * phi for phi in w ], c

m = 10
n = 9
u, c = u_fun(x, y, m, n)



## Matrices and solution



In [8]:

dudx = [diff(u[k], x) for k in range(len(c))]
dudy = [diff(u[k], y) for k in range(len(c))]




In [9]:

Kaux = Matrix(m*n, m*n, lambda ii, jj: dudx[ii]*dudx[jj] + dudy[ii]*dudy[jj])
Maux = Matrix(m*n, m*n, lambda ii, jj: u[ii]*u[jj])
K = Matrix(m*n, m*n, lambda i,j: 0)
M = Matrix(m*n, m*n, lambda i,j: 0)



The integrals should be of the form

$$B_{ij} = \int\limits_{-1}^1\int\limits_{\frac{1}{2}(x^2 - 1)}^{\frac{1}{2}(1 - x^2)} A_{ij} \mathrm{d}y\, \mathrm{d}x\, ,$$


In [10]:

for row in range(m*n):
for col in range(row, m*n):
K_inte = Kaux[row, col]
M_inte = Maux[row, col]
K_inte = integrate(K_inte, (y, (x**2 - 1)/2, (1 - x**2)/2), (x, -1, 1))
M_inte = integrate(M_inte, (y, (x**2 - 1)/2, (1 - x**2)/2), (x, -1, 1))
K[row, col] += K_inte
M[row, col] += M_inte
if row != col:
K[col, row] += K_inte
M[col, row] += M_inte




In [11]:

Kn = np.array(K).astype(np.float64)
Mn = np.array(M).astype(np.float64)




In [12]:

vals, vecs = eigh(Kn, Mn, eigvals=(0,8))
vals




Out[12]:

array([ 16.10095349,  30.93334613,  50.33837571,  50.33838703,
74.38718728,  74.38735248, 102.69815051, 102.75135697,
105.00010187])



## Visualization of the modes



In [13]:




In [14]:

plt.figure(figsize=(8, 8))
for i in range(8):
U = sum(vecs[j, i]*u[j] for j in range(m*n))
vecU = lambdify((x,y), U, "numpy")
Z_max = Z.max()
Z_max = max (Z_max, -Z.min())
plt.subplot(4, 2, i + 1)
plt.title(r"$k^2=%.2f$" % vals[i], size=12);
plt.contour(X, Y, bound_num, [0], colors="black")
plt.contourf(X, Y, Z, 12, cmap="RdYlBu", vmin=-1.2, vmax=1.2)
plt.axis("image")
plt.axis(False)
plt.savefig("../../../Documents/%dx%d.png" % (m, n), dpi=300)




var element = $('#ca1ea389-3348-4035-851b-9f36496b658a'); /* Put everything inside the global mpl namespace */ window.mpl = {}; mpl.get_websocket_type = function() { if (typeof(WebSocket) !== 'undefined') { return WebSocket; } else if (typeof(MozWebSocket) !== 'undefined') { return MozWebSocket; } else { alert('Your browser does not have WebSocket support. ' + 'Please try Chrome, Safari or Firefox ≥ 6. ' + 'Firefox 4 and 5 are also supported but you ' + 'have to enable WebSockets in about:config.'); }; } mpl.figure = function(figure_id, websocket, ondownload, parent_element) { this.id = figure_id; this.ws = websocket; this.supports_binary = (this.ws.binaryType != undefined); if (!this.supports_binary) { var warnings = document.getElementById("mpl-warnings"); if (warnings) { warnings.style.display = 'block'; warnings.textContent = ( "This browser does not support binary websocket messages. " + "Performance may be slow."); } } this.imageObj = new Image(); this.context = undefined; this.message = undefined; this.canvas = undefined; this.rubberband_canvas = undefined; this.rubberband_context = undefined; this.format_dropdown = undefined; this.image_mode = 'full'; this.root =$('<div/>');
this._root_extra_style(this.root)
this.root.attr('style', 'display: inline-block');

$(parent_element).append(this.root); this._init_header(this); this._init_canvas(this); this._init_toolbar(this); var fig = this; this.waiting = false; this.ws.onopen = function () { fig.send_message("supports_binary", {value: fig.supports_binary}); fig.send_message("send_image_mode", {}); if (mpl.ratio != 1) { fig.send_message("set_dpi_ratio", {'dpi_ratio': mpl.ratio}); } fig.send_message("refresh", {}); } this.imageObj.onload = function() { if (fig.image_mode == 'full') { // Full images could contain transparency (where diff images // almost always do), so we need to clear the canvas so that // there is no ghosting. fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height); } fig.context.drawImage(fig.imageObj, 0, 0); }; this.imageObj.onunload = function() { fig.ws.close(); } this.ws.onmessage = this._make_on_message_function(this); this.ondownload = ondownload; } mpl.figure.prototype._init_header = function() { var titlebar =$(
'<div class="ui-dialog-titlebar ui-widget-header ui-corner-all ' +
'ui-helper-clearfix"/>');
var titletext = ( '<div class="ui-dialog-title" style="width: 100%; ' + 'text-align: center; padding: 3px;"/>'); titlebar.append(titletext) this.root.append(titlebar); this.header = titletext[0]; } mpl.figure.prototype._canvas_extra_style = function(canvas_div) { } mpl.figure.prototype._root_extra_style = function(canvas_div) { } mpl.figure.prototype._init_canvas = function() { var fig = this; var canvas_div =('<div/>');

canvas_div.attr('style', 'position: relative; clear: both; outline: 0');

function canvas_keyboard_event(event) {
return fig.key_event(event, event['data']);
}

canvas_div.keydown('key_press', canvas_keyboard_event);
canvas_div.keyup('key_release', canvas_keyboard_event);
this.canvas_div = canvas_div
this._canvas_extra_style(canvas_div)
this.root.append(canvas_div);

var canvas = $('<canvas/>'); canvas.addClass('mpl-canvas'); canvas.attr('style', "left: 0; top: 0; z-index: 0; outline: 0") this.canvas = canvas[0]; this.context = canvas[0].getContext("2d"); var backingStore = this.context.backingStorePixelRatio || this.context.webkitBackingStorePixelRatio || this.context.mozBackingStorePixelRatio || this.context.msBackingStorePixelRatio || this.context.oBackingStorePixelRatio || this.context.backingStorePixelRatio || 1; mpl.ratio = (window.devicePixelRatio || 1) / backingStore; var rubberband =$('<canvas/>');
rubberband.attr('style', "position: absolute; left: 0; top: 0; z-index: 1;")

var pass_mouse_events = true;

canvas_div.resizable({
start: function(event, ui) {
pass_mouse_events = false;
},
resize: function(event, ui) {
fig.request_resize(ui.size.width, ui.size.height);
},
stop: function(event, ui) {
pass_mouse_events = true;
fig.request_resize(ui.size.width, ui.size.height);
},
});

function mouse_event_fn(event) {
if (pass_mouse_events)
return fig.mouse_event(event, event['data']);
}

rubberband.mousedown('button_press', mouse_event_fn);
rubberband.mouseup('button_release', mouse_event_fn);
// Throttle sequential mouse events to 1 every 20ms.
rubberband.mousemove('motion_notify', mouse_event_fn);

rubberband.mouseenter('figure_enter', mouse_event_fn);
rubberband.mouseleave('figure_leave', mouse_event_fn);

canvas_div.on("wheel", function (event) {
event = event.originalEvent;
event['data'] = 'scroll'
if (event.deltaY < 0) {
event.step = 1;
} else {
event.step = -1;
}
mouse_event_fn(event);
});

canvas_div.append(canvas);
canvas_div.append(rubberband);

this.rubberband = rubberband;
this.rubberband_canvas = rubberband[0];
this.rubberband_context = rubberband[0].getContext("2d");
this.rubberband_context.strokeStyle = "#000000";

this._resize_canvas = function(width, height) {
// Keep the size of the canvas, canvas container, and rubber band
// canvas in synch.
canvas_div.css('width', width)
canvas_div.css('height', height)

canvas.attr('width', width * mpl.ratio);
canvas.attr('height', height * mpl.ratio);
canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');

rubberband.attr('width', width);
rubberband.attr('height', height);
}

// Set the figure to an initial 600x600px, this will subsequently be updated
// upon first draw.
this._resize_canvas(600, 600);

// Disable right mouse context menu.
$(this.rubberband_canvas).bind("contextmenu",function(e){ return false; }); function set_focus () { canvas.focus(); canvas_div.focus(); } window.setTimeout(set_focus, 100); } mpl.figure.prototype._init_toolbar = function() { var fig = this; var nav_element =$('<div/>');
nav_element.attr('style', 'width: 100%');
this.root.append(nav_element);

// Define a callback function for later on.
function toolbar_event(event) {
return fig.toolbar_button_onclick(event['data']);
}
function toolbar_mouse_event(event) {
return fig.toolbar_button_onmouseover(event['data']);
}

for(var toolbar_ind in mpl.toolbar_items) {
var name = mpl.toolbar_items[toolbar_ind][0];
var tooltip = mpl.toolbar_items[toolbar_ind][1];
var image = mpl.toolbar_items[toolbar_ind][2];
var method_name = mpl.toolbar_items[toolbar_ind][3];

if (!name) {
// put a spacer in here.
continue;
}
var button = $('<button/>'); button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' + 'ui-button-icon-only'); button.attr('role', 'button'); button.attr('aria-disabled', 'false'); button.click(method_name, toolbar_event); button.mouseover(tooltip, toolbar_mouse_event); var icon_img =$('<span/>');

var tooltip_span = $('<span/>'); tooltip_span.addClass('ui-button-text'); tooltip_span.html(tooltip); button.append(icon_img); button.append(tooltip_span); nav_element.append(button); } var fmt_picker_span =$('<span/>');

var fmt_picker = $('<select/>'); fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content'); fmt_picker_span.append(fmt_picker); nav_element.append(fmt_picker_span); this.format_dropdown = fmt_picker[0]; for (var ind in mpl.extensions) { var fmt = mpl.extensions[ind]; var option =$(
'<option/>', {selected: fmt === mpl.default_extension}).html(fmt);
fmt_picker.append(option);
}

// Add hover states to the ui-buttons
$( ".ui-button" ).hover( function() {$(this).addClass("ui-state-hover");},
function() { $(this).removeClass("ui-state-hover");} ); var status_bar =$('<span class="mpl-message"/>');
nav_element.append(status_bar);
this.message = status_bar[0];
}

mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {
// Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,
// which will in turn request a refresh of the image.
this.send_message('resize', {'width': x_pixels, 'height': y_pixels});
}

mpl.figure.prototype.send_message = function(type, properties) {
properties['type'] = type;
properties['figure_id'] = this.id;
this.ws.send(JSON.stringify(properties));
}

mpl.figure.prototype.send_draw_message = function() {
if (!this.waiting) {
this.waiting = true;
this.ws.send(JSON.stringify({type: "draw", figure_id: this.id}));
}
}

mpl.figure.prototype.handle_save = function(fig, msg) {
var format_dropdown = fig.format_dropdown;
var format = format_dropdown.options[format_dropdown.selectedIndex].value;
}

mpl.figure.prototype.handle_resize = function(fig, msg) {
var size = msg['size'];
if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {
fig._resize_canvas(size[0], size[1]);
fig.send_message("refresh", {});
};
}

mpl.figure.prototype.handle_rubberband = function(fig, msg) {
var x0 = msg['x0'] / mpl.ratio;
var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;
var x1 = msg['x1'] / mpl.ratio;
var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;
x0 = Math.floor(x0) + 0.5;
y0 = Math.floor(y0) + 0.5;
x1 = Math.floor(x1) + 0.5;
y1 = Math.floor(y1) + 0.5;
var min_x = Math.min(x0, x1);
var min_y = Math.min(y0, y1);
var width = Math.abs(x1 - x0);
var height = Math.abs(y1 - y0);

fig.rubberband_context.clearRect(
0, 0, fig.canvas.width / mpl.ratio, fig.canvas.height / mpl.ratio);

fig.rubberband_context.strokeRect(min_x, min_y, width, height);
}

mpl.figure.prototype.handle_figure_label = function(fig, msg) {
}

mpl.figure.prototype.handle_cursor = function(fig, msg) {
var cursor = msg['cursor'];
switch(cursor)
{
case 0:
cursor = 'pointer';
break;
case 1:
cursor = 'default';
break;
case 2:
cursor = 'crosshair';
break;
case 3:
cursor = 'move';
break;
}
fig.rubberband_canvas.style.cursor = cursor;
}

mpl.figure.prototype.handle_message = function(fig, msg) {
fig.message.textContent = msg['message'];
}

mpl.figure.prototype.handle_draw = function(fig, msg) {
// Request the server to send over a new figure.
fig.send_draw_message();
}

mpl.figure.prototype.handle_image_mode = function(fig, msg) {
fig.image_mode = msg['mode'];
}

mpl.figure.prototype.updated_canvas_event = function() {
// Called whenever the canvas gets updated.
this.send_message("ack", {});
}

// A function to construct a web socket function for onmessage handling.
// Called in the figure constructor.
mpl.figure.prototype._make_on_message_function = function(fig) {
return function socket_on_message(evt) {
if (evt.data instanceof Blob) {
/* FIXME: We get "Resource interpreted as Image but
* transferred with MIME type text/plain:" errors on
* Chrome.  But how to set the MIME type?  It doesn't seem
* to be part of the websocket stream */
evt.data.type = "image/png";

/* Free the memory for the previous frames */
if (fig.imageObj.src) {
(window.URL || window.webkitURL).revokeObjectURL(
fig.imageObj.src);
}

fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(
evt.data);
fig.updated_canvas_event();
fig.waiting = false;
return;
}
else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == "data:image/png;base64") {
fig.imageObj.src = evt.data;
fig.updated_canvas_event();
fig.waiting = false;
return;
}

var msg = JSON.parse(evt.data);
var msg_type = msg['type'];

// Call the  "handle_{type}" callback, which takes
// the figure and JSON message as its only arguments.
try {
var callback = fig["handle_" + msg_type];
} catch (e) {
console.log("No handler for the '" + msg_type + "' message type: ", msg);
return;
}

if (callback) {
try {
// console.log("Handling '" + msg_type + "' message: ", msg);
callback(fig, msg);
} catch (e) {
console.log("Exception inside the 'handler_" + msg_type + "' callback:", e, e.stack, msg);
}
}
};
}

// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas
mpl.findpos = function(e) {
//this section is from http://www.quirksmode.org/js/events_properties.html
var targ;
if (!e)
e = window.event;
if (e.target)
targ = e.target;
else if (e.srcElement)
targ = e.srcElement;
if (targ.nodeType == 3) // defeat Safari bug
targ = targ.parentNode;

// jQuery normalizes the pageX and pageY
// pageX,Y are the mouse positions relative to the document
// offset() returns the position of the element relative to the document
var x = e.pageX - $(targ).offset().left; var y = e.pageY -$(targ).offset().top;

return {"x": x, "y": y};
};

/*
* return a copy of an object with only non-object keys
* we need this to avoid circular references
* http://stackoverflow.com/a/24161582/3208463
*/
function simpleKeys (original) {
return Object.keys(original).reduce(function (obj, key) {
if (typeof original[key] !== 'object')
obj[key] = original[key]
return obj;
}, {});
}

mpl.figure.prototype.mouse_event = function(event, name) {
var canvas_pos = mpl.findpos(event)

if (name === 'button_press')
{
this.canvas.focus();
this.canvas_div.focus();
}

var x = canvas_pos.x * mpl.ratio;
var y = canvas_pos.y * mpl.ratio;

this.send_message(name, {x: x, y: y, button: event.button,
step: event.step,
guiEvent: simpleKeys(event)});

/* This prevents the web browser from automatically changing to
* the text insertion cursor when the button is pressed.  We want
* to control all of the cursor setting manually through the
* 'cursor' event from matplotlib */
event.preventDefault();
return false;
}

mpl.figure.prototype._key_event_extra = function(event, name) {
// Handle any extra behaviour associated with a key event
}

mpl.figure.prototype.key_event = function(event, name) {

// Prevent repeat events
if (name == 'key_press')
{
if (event.which === this._key)
return;
else
this._key = event.which;
}
if (name == 'key_release')
this._key = null;

var value = '';
if (event.ctrlKey && event.which != 17)
value += "ctrl+";
if (event.altKey && event.which != 18)
value += "alt+";
if (event.shiftKey && event.which != 16)
value += "shift+";

value += 'k';
value += event.which.toString();

this._key_event_extra(event, name);

this.send_message(name, {key: value,
guiEvent: simpleKeys(event)});
return false;
}

mpl.figure.prototype.toolbar_button_onclick = function(name) {
this.handle_save(this, null);
} else {
this.send_message("toolbar_button", {name: name});
}
};

mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {
this.message.textContent = tooltip;
};
mpl.toolbar_items = [["Home", "Reset original view", "fa fa-home icon-home", "home"], ["Back", "Back to previous view", "fa fa-arrow-left icon-arrow-left", "back"], ["Forward", "Forward to next view", "fa fa-arrow-right icon-arrow-right", "forward"], ["", "", "", ""], ["Pan", "Pan axes with left mouse, zoom with right", "fa fa-arrows icon-move", "pan"], ["Zoom", "Zoom to rectangle", "fa fa-square-o icon-check-empty", "zoom"], ["", "", "", ""], ["Download", "Download plot", "fa fa-floppy-o icon-save", "download"]];

mpl.extensions = ["eps", "jpeg", "pdf", "png", "ps", "raw", "svg", "tif"];

mpl.default_extension = "png";var comm_websocket_adapter = function(comm) {
// Create a "websocket"-like object which calls the given IPython comm
// object with the appropriate methods. Currently this is a non binary
// socket, so there is still some room for performance tuning.
var ws = {};

ws.close = function() {
comm.close()
};
ws.send = function(m) {
//console.log('sending', m);
comm.send(m);
};
// Register the callback with on_msg.
comm.on_msg(function(msg) {
//console.log('receiving', msg['content']['data'], msg);
// Pass the mpl event to the overridden (by mpl) onmessage function.
ws.onmessage(msg['content']['data'])
});
return ws;
}

mpl.mpl_figure_comm = function(comm, msg) {
// This is the function which gets called when the mpl process
// starts-up an IPython Comm through the "matplotlib" channel.

var id = msg.content.data.id;
// Get hold of the div created by the display call when the Comm
// socket was opened in Python.
var element = $("#" + id); var ws_proxy = comm_websocket_adapter(comm) function ondownload(figure, format) { window.open(figure.imageObj.src); } var fig = new mpl.figure(id, ws_proxy, ondownload, element.get(0)); // Call onopen now - mpl needs it, as it is assuming we've passed it a real // web socket which is closed, not our websocket->open comm proxy. ws_proxy.onopen(); fig.parent_element = element.get(0); fig.cell_info = mpl.find_output_cell("<div id='" + id + "'></div>"); if (!fig.cell_info) { console.error("Failed to find cell for figure", id, fig); return; } var output_index = fig.cell_info[2] var cell = fig.cell_info[0]; }; mpl.figure.prototype.handle_close = function(fig, msg) { var width = fig.canvas.width/mpl.ratio fig.root.unbind('remove') // Update the output cell to use the data from the current canvas. fig.push_to_output(); var dataURL = fig.canvas.toDataURL(); // Re-enable the keyboard manager in IPython - without this line, in FF, // the notebook keyboard shortcuts fail. IPython.keyboard_manager.enable()$(fig.parent_element).html('<amp-img layout="responsive" width="500" height="300" src="' + dataURL + '" width="' + width + '">');
fig.close_ws(fig, msg);
}

mpl.figure.prototype.close_ws = function(fig, msg){
fig.send_message('closing', msg);
// fig.ws.close()
}

mpl.figure.prototype.push_to_output = function(remove_interactive) {
// Turn the data on the canvas into data in the output cell.
var width = this.canvas.width/mpl.ratio
var dataURL = this.canvas.toDataURL();
this.cell_info[1]['text/html'] = '<amp-img layout="responsive" width="500" height="300" src="' + dataURL + '" width="' + width + '">';
}

mpl.figure.prototype.updated_canvas_event = function() {
// Tell IPython that the notebook contents must change.
IPython.notebook.set_dirty(true);
this.send_message("ack", {});
var fig = this;
// Wait a second, then push the new image to the DOM so
// that it is saved nicely (might be nice to debounce this).
setTimeout(function () { fig.push_to_output() }, 1000);
}

mpl.figure.prototype._init_toolbar = function() {
var fig = this;

var nav_element = $('<div/>'); nav_element.attr('style', 'width: 100%'); this.root.append(nav_element); // Define a callback function for later on. function toolbar_event(event) { return fig.toolbar_button_onclick(event['data']); } function toolbar_mouse_event(event) { return fig.toolbar_button_onmouseover(event['data']); } for(var toolbar_ind in mpl.toolbar_items){ var name = mpl.toolbar_items[toolbar_ind][0]; var tooltip = mpl.toolbar_items[toolbar_ind][1]; var image = mpl.toolbar_items[toolbar_ind][2]; var method_name = mpl.toolbar_items[toolbar_ind][3]; if (!name) { continue; }; var button =$('<button class="btn btn-default" href="#" title="' + name + '"><i class="fa ' + image + ' fa-lg"></i></button>');
button.click(method_name, toolbar_event);
button.mouseover(tooltip, toolbar_mouse_event);
nav_element.append(button);
}

var status_bar = ('<span class="mpl-message" style="text-align:right; float: right;"/>'); nav_element.append(status_bar); this.message = status_bar[0]; // Add the close button to the window. var buttongrp =('<div class="btn-group inline pull-right"></div>');
var button = $('<button class="btn btn-mini btn-primary" href="#" title="Stop Interaction"><i class="fa fa-power-off icon-remove icon-large"></i></button>'); button.click(function (evt) { fig.handle_close(fig, {}); } ); button.mouseover('Stop Interaction', toolbar_mouse_event); buttongrp.append(button); var titlebar = this.root.find($('.ui-dialog-titlebar'));
titlebar.prepend(buttongrp);
}

mpl.figure.prototype._root_extra_style = function(el){
var fig = this
el.on("remove", function(){
fig.close_ws(fig, {});
});
}

mpl.figure.prototype._canvas_extra_style = function(el){
// this is important to make the div 'focusable
el.attr('tabindex', 0)
// reach out to IPython and tell the keyboard manager to turn it's self
// off when our div gets focus

// location in version 3
if (IPython.notebook.keyboard_manager) {
IPython.notebook.keyboard_manager.register_events(el);
}
else {
// location in version 2
IPython.keyboard_manager.register_events(el);
}

}

mpl.figure.prototype._key_event_extra = function(event, name) {
var manager = IPython.notebook.keyboard_manager;
if (!manager)
manager = IPython.keyboard_manager;

// Check for shift+enter
if (event.shiftKey && event.which == 13) {
this.canvas_div.blur();
// select the cell after this one
var index = IPython.notebook.find_cell_index(this.cell_info[0]);
IPython.notebook.select(index + 1);
}
}

mpl.figure.prototype.handle_save = function(fig, msg) {
}

mpl.find_output_cell = function(html_output) {
// Return the cell and output element which can be found *uniquely* in the notebook.
// Note - this is a bit hacky, but it is done because the "notebook_saving.Notebook"
// IPython event is triggered only after the cells have been serialised, which for
// our purposes (turning an active figure into a static one), is too late.
var cells = IPython.notebook.get_cells();
var ncells = cells.length;
for (var i=0; i<ncells; i++) {
var cell = cells[i];
if (cell.cell_type === 'code'){
for (var j=0; j<cell.output_area.outputs.length; j++) {
var data = cell.output_area.outputs[j];
if (data.data) {
// IPython >= 3 moved mimebundle to data attribute of output
data = data.data;
}
if (data['text/html'] == html_output) {
return [cell, data, j];
}
}
}
}
}

// Register the function which deals with the matplotlib target/channel.
// The kernel may be null if the page has been refreshed.
if (IPython.notebook.kernel != null) {
IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);
}



## Analytic mass matrix

The mass matrix can be integrated analytically.



In [15]:

from sympy import gamma

def mass_coeff(j, k, m, n):
coeff = (k + n + 2)*(k + n + 4)*(1 + (-1)**(j + m))*(1 + (-1)**(k + n))
coeff *= gamma(n + k + 1)*gamma((1 + m + j)/2)/gamma((13 + m + j + 2*n + 2*k)/2)
coeff /= 2**(k + n - 1)
return coeff




In [16]:

mass_coeff(0, 0, 0, 0)




Out[16]:

$\displaystyle 0.394035594035594$




In [20]:

vals




Out[20]:

array([ 16.10095349,  30.93334613,  50.33837571,  50.33838703,
74.38718728,  74.38735248, 102.69815051, 102.75135697,
105.00010187])




In [ ]: