Matlab: Resolviendo la ecuación de onda bidimensional en el caso de una membrana de forma piramidal

El código fuente para resolver una ecuación de onda bidimensional está dado en el archivo twodimwavedirbc.m disponible en esta web: MATH 5343 Numerical Solution for Partial Differential Equations, Spring 2011.

La idea es resolver el ejercicio 12a del capítulo 12.2 del libro "Introduction to Numerical Ordinary and Partial Differential Equations Using MATLAB" de Alexander Stanoyevitch (de este libro proviene el scrip twodimwavedirbc.m):

Es decir, es una ecuación de onda cuya forma inicial es una pirámide. La dificultad aquí es dibujar la pirámide, yo la deduje como explico a continuacion:

La pirámide mostrada:

está formada por cuatro planos limitados por las rectas Y=X y Y = -X/2+4:

Las regiones donde Z es una función de X o de Y (dependiendo del plano donde se proyecta el eje Z) se muestran a continuación:

Para X<2, Y<2 y la región donde X<Y, la variable Z está definida por la función Z= Y/2.
Esto se hace para cada región de cada uno de los cuatro planos que definen la pirámide.

El script resultante es:

function z = pyramid(x, y)
z=0;

if y<=2 && x <= 2
   if y < x
       z=y/2;
   else
       z=x/2;
   end
end

if y<=2 && x > 2 && x <= 4
   if y < 4-x
       z=y/2;
   else
       z=-x/2+2;
   end
end

if y > 2 && x <= 2 && y <= 4
   if y < 4-x
       z=x/2;
   else
       z=-y/2+2;
   end
end

if y > 2 && y <= 4 && x > 2 && x <= 4
   if y < x
       z=-x/2+2;
   else
       z=-y/2+2;
   end
end

end 

Modifiqué un poco el archivo twodimwavedirbc.m para hacerlo más flexible:

function [x, y, t, U] = twodimwavedirbc2(phi, nu, a1, a, b1, b, T, h, c)
% a1: punto inicial del eje x
% b1: punto inicial del eje y

% solves the two-dimensional wave problem u_tt = c^2(u_xx+u_yy)
% on the rectangle {a1<= x <= a,  b1<= y <= b}, with u(x,y)=0
% on the boundary of the rectangle.

% variables de entrada: 
% phi=phi(x,y) = función de la forma inicial de la onda
% nu=nu(x,y) = función inicial de la velocidad onda

% a= right endpoint of x, b = right endpoint of y, 
% T = final time solution will be computed.  h=common gap on
% x, y-grids,   k= t-grid gap, c = speed of wave.
% wave.
% Output variables: x = row vector for first space variable, y =
% row vector for second space variable,  t = time grid row vector
%(starts at t=0, ends at t=T, has Nt equally spaced values), 
% U = (Nx)by(Ny)by(Nt) matrix of solution approximations at
% corresponding grid points (where Ny = number of y-grid points)
% x grid will correspond to first entries of U, y
% grid values to second entries of U, and t grid to third entries of U.
% CAUTION:  This method will only work if h is chosen so that the x
% and y grids can have a common gap size, i.e., if h = a/h and
% b/h must be integers.
% The time grid gap is chosen so that mu = 1/2, this guarantees the
% Courant-Friedrichs-Levy condition holds and simplifies the
% main finite difference formula.

k=h/sqrt(2)/c; %k is determined from mu = 1/2

if (b/h-floor(round(b/h))>10*eps)||(a/h-floor(round(a/h))>10*eps)
    fprintf('Space grid gap h must divide evenly into both a and b \r')
    fprintf(' Either try another input or modify the algorithm')
    error('M-file will exit')
end

Nx = round((a-a1)/h)+1; %number of points on x-grid
Ny = round((b-b1)/h)+1; %number of points on y-grid
Nt = floor(T/k)+1; %number of points on t-grid
U=zeros(Nx, Ny, Nt);
x=a1:h:a;
y=b1:h:b;
t=0:k:T;
% Recall matrix indices must start at 1.  Thus the indices of the
% matrix will always be one more than the corresponding indices that
% were used in theoretical development.

% Note that by default, zero boundary values have been assigned to all
% grid points on the edges of the rectangle (and, for the time being,
% at all grid points).
%Assign initial time t=0 values and next step t=k values.
for i=2:(Nx-1)
     for j=2:(Ny-1)
    U(i,j,1)=feval(phi,x(i),y(j));
    U(i,j,2)=.25*(feval(phi,x(i-1),y(j))+...
      feval(phi,x(i+1),y(j))+feval(phi,x(i),y(j-1))+...
      feval(phi,x(i),y(j+1)))+ k*feval(nu,x(i),y(j));     
     end
end
%Assign values at interior grid points
for ell=3:Nt %letter ell looks too much like number one
    U(2:(Nx-1),2:(Ny-1), ell) = ...
     +.5*(U(3:Nx,2:(Ny-1), ell-1)+ U(1:(Nx-2),2:(Ny-1), ell-1)...
     + U(2:(Nx-1),3:Ny, ell-1) + U(2:(Nx-1),1:(Ny-2), ell-1))...
     - U(2:(Nx-1),2:(Ny-1), ell-2);
end

Y el script principal es:

clc; clear all; close all;
phi = @(x, y)sin(2*x) *sin(2*y);
nu=@(x, y)0;
[x, y, t, U] = twodimwavedirbc2('pyramid', nu, 0, 4, 0, 4, 4, 0.1, 1);
[a b c] = size(U);
%ans = 26 26 46
for ell=1:c
     surf(x,y,U(:, :,ell));
     axis([0, 4, 0, 4 -1, 1]);
     M(:,:,ell) = getframe;
end
movie(M, 2, 4)


Lo he probado en Matlab, y acá está corriendo lo mejor que puede en Octave 4.4.1 a pesar de que me dice que la función Movie no está implementada:

Resolviendo una ecuación parabólica unidimensional en Matlab/Octave

El siguiente ejercicio se encuentra en el libro "Computational Partial Differential Equations Using MATLAB" by Yi-Tung Chen, Jichun Li, capítulo 2, ejercicio 3:

Una pista de cómo resolver este ejercicio está en la sección 2.21:

Los autores dan un código base sobre el qué trabajar. Luego de recibir un poco de ayuda de la gentita de Stackoverflow, este es el código para resolver el problema:


 %---------------------------------------------------------------
% para1d.m:
%    use the explicit scheme to solve the parabolic equation
%    u_t(x,t) = u_{xx}(x,t),          xl < x < xr, 0 < t < tf
%    u(x,0) = f(x),                   xl < x < xr
%    u_x(0,t) = gl(t), u_x(1,t) = gr(t),  0  < t < tf
%---------------------------------------------------------------
clc;
close all;
clear all;                  % clear all variables in memory

xl=0; xr=1;                 % x domain [xl,xr]
J = 31;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 51;                    % Nt: number of time steps
dt = tf/Nt/4;

mu = dt/(dx)^2;

if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end

% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
f = x;

% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
u_ex = zeros(J+1,Nt);

% boundary condition at left side
gl = 0;
% boundary condition at right side
gr = 0;

% Find the approximate solution at each time step
for n = 1:Nt
    t = n*dt;         % current time
     
    if n==1
      for j=2:J
       % first time step
        u(j,1) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
      end
    end
  
    if n > 1
        for j=2:J
         % interior nodes
           u(j,n) = u(j,n-1) + mu*(u(j+1,n-1) - 2*u(j,n-1) + u(j-1,n-1));
        end
    end
  
    % hint 2.21 Los valores en la frontera es mejor calcularlos al final
    % pues estoy cogiendo índices que dependen de los nodos interiores
    u(1, n) = u(2, n) - dx * gl;
    u(J+1,n) = u(J, n) + dx * gr;
  
    % calculate the analytic solution
    for j=1:J+1
      xj = xl + (j-1)*dx;    
      suma = zeros(1000 , 1);

      for k= 1:1000
        suma(k) = 4/(((2*k-1)*pi)^2);
        suma(k) = suma(k) * exp(-((2*k-1)^2)*t*pi^2) * cos((2*k-1)*pi*xj);
      end

      u_ex(j, n)= double(0.5 - sum(suma));
    end
end

% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(cool);     % draw gray figure
surf(x, tt, u');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')

figure(2)
surf(x, tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')

maxerr=max(max(abs(u-u_ex))),


El resultado es:

Probado en Octave 4.4.1 y Matlab 2008 R2.