rsvp_ss18_full

rsvp_ss18_full.pdf PDF document, 1.65 MB

heat_mol.py

heat_mol.py Python Source, 2 KB

File contents

#----------------------------------------------------------------------
#Python example code for solving the heat equation with boundary input
#using method of lines: finite differences with dopri5
#winter term 2020
#Alexander Schaum
#----------------
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import ode
import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 22})
import numpy as np
import csv

plt.rc('text', usetex=True)
plt.rc('font', family='serif')

#number of discretization points in space, space vector, discretization step
N  = 40
z = np.linspace(0,1,N)
dz = z[1]-z[0]
#simulation time, number of time points to be included in solution vector, solution time vector, time step
T = 5;
tN = 400
t = np.linspace(0,T,tN)
dt = t[1]-t[0]
#coefficient
a = 0.5
#save figure option (YES=1,NO=0)
save = 0

def rhs(t, x):
  #decaying periodic exitation signal
  v = np.sin(5.0*t)*np.exp(-1.0*t)
  return A.dot(x) + f(x) + b.dot(v)

def f(x):
  #returns the source term
  return 3.0*x*(1.0-x)

#matrices for numerical computation (discretized model equations)  
A = (-2*np.eye(N) + np.diag(np.ones(N-1),k=-1) + np.diag(np.ones(N-1),k=1))
A[0,0] = -1
A = a/(dz**2.0)*A

b = np.zeros(N)
b[N-1]=a/(dz**2.0)
#setting up the numerical ODE integration method (dopri5)
x0 = np.cos(z*np.pi/2.0)
ss = ode(lambda t,x: rhs(t,x), jac=None)
ss.set_integrator('dopri5',method='bdf')
ss.set_f_params()
ss.set_initial_value(x0, 0.0)
#defining solution matrix (space,time)
sol = np.zeros(shape=(N,tN))
sol[:,0]=x0
#solving the system of ODEs
for i in range(1,tN):
  ss.integrate(ss.t+dt)
  sol[:,i] = ss.y
#graphical output
T, Z = np.meshgrid(t,z)
fig = plt.figure(figsize=(10,6))
fig.set_facecolor('white')
ax = fig.add_subplot(1,1,1,projection='3d')
ax.w_xaxis.set_pane_color((1.0,1.0,1.0,1.0))
ax.w_yaxis.set_pane_color((1.0,1.0,1.0,1.0))
ax.w_zaxis.set_pane_color((1.0,1.0,1.0,1.0))
p = ax.plot_surface(T,Z,sol,cmap=cm.plasma)
ax.set_xlabel(r'$t$', fontsize=14)
ax.set_ylabel(r'$z$', fontsize=14, rotation=0)
ax.set_zlabel(r'$x(z,t)$', fontsize=14)
ax.view_init(37,111)#elevation, azimuth
#optional: save the figure (e.g., as .eps)
if save == 1:
  plt.savefig("heat_mol.eps",bbox_inches='tight', pad_inches = 0)

plt.show()

ex_pdepe

ex_pdepe.m Objective-C source code, 1 KB

File contents

function ex1_dr_sys()
%Example of pdepe function usage for a diffusion-convection-reaction system
%----------------------------------
% pdt x = a pdz^2 x - v pdz x + k x
% x(0,t)=0, pdz x(1,t)=u
% x(z,0) = sin(pi/2 z)
%-----------------------
%in pdepe notation:
% c*pdt x = x^(-m)\pdz f(x,DxDz) + s(x,DxDz)
% left boudnary
% ql+pl*f -ul=0
% right boundary
% qr+pr*f -ur=0
% yielding: c=1, f=a DxDz,, s=-v*Dxdz + kx
%           pl=xl, ql=0
%           pr=-u, qr=1/a
%(c 2020) Alexander Schaum, CAU Kiel
%-------------------------------------------------------------------- 
close all; 

%--------------------------------------------------------- 
%spatial coordinate discretization points
z = 0:0.01:1.0; 

%system parameters
p.a = 0.5;  %diffusion coefficient
p.v = 0.5;  %velocity of advective flow
p.k = 1.75;  %reaction rate

%---------------------------------------------------------
%Simulation using pdepe
m = 0;
z = linspace(z(1),z(end),101);
t=linspace(0,20,200);

%Solver call
sol = pdepe(m,@dcrs_pde,@dcrs_ic,@dcrs_bc,z,t,[],p);
sol = sol(:,:,1);

figure(112);
mesh(z,t,sol);
az = -70;
el = 40;
view(az,el)
xlabel('z')
ylabel('t')
zlabel('x(z,t)')

%---------------------------------------------------------
% SUBFUNCTIONS

function [c,f,s]=dcrs_pde(z,t,x,DxDz,p)
%Partial differential equation
c = 1.0;
f = p.a*DxDz;
s = -p.v*DxDz + p.k*x;

function x0=dcrs_ic(z,p);
%Initial condition
x0 = sin(pi/2*z); 

function [pl,ql,pr,qr]=dcrs_bc(zl,xl,zr,xr,t,p)
%Boundary conditions
u  = sin(2*t)*exp(-0.1*t);
pl = xl;
ql = 0;
pr = -u;
qr = 1/p.a;

heat_mol_matlab

heat_mol.m Objective-C source code, 1 KB

File contents

function heat_mol()

a = 0.5;%"diffusion" coefficient

N=40;%number of discretization steps
z = linspace(0,1,N);%discretized space interval
dz = z(2)-z(1);%spatial step size
T = 5;%simultion time
tN = 200;%number of interpolation points
t = linspace(0,T,tN);%time vector for numerical solution output
%system matrixes
A = (-2*eye(N) + diag(ones(1,N-1),-1) + diag(ones(1,N-1),1));
A(1,1) = -1;
p.A = a/(dz^2.0)*A;

b = zeros(N,1);
b(N)=a/(dz^2.0);
p.b = b;

%initial condition and numerical solution
x0 = sin(pi/2*z).^2;
[tn,xn]=ode45(@(t,x)rhs(t,x,p),t,x0);
%graphical output
figure(1)
mesh(z,tn,xn)
az = -70;
el = 40;
view(az,el)
xlabel('z')
ylabel('t')
zlabel('x(z,t)')

end

function out = rhs(t,x,p)
 v = sin(5.0*t)*exp(-1.0*t);
 out = p.A*x+p.b*v;
end
Research