# hyperbolic

Solve hyperbolic PDE problem

Hyperbolic equation solver

Solves PDE problems of the type

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

on a 2-D or 3-D region Ω, or the system PDE problem

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

The variables c, a, f, and d can depend on position, time, and the solution u and its gradient.

## Syntax

• u = hyperbolic(u0,ut0,tlist,model,c,a,f,d)
example
• u = hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)
example
• u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M)
example
• u = hyperbolic(___,rtol)
• u = hyperbolic(___,rtol,atol)
• u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M,___,'DampingMatrix',D)
example
• u = hyperbolic(___,'Stats','off')

## Description

example

u = hyperbolic(u0,ut0,tlist,model,c,a,f,d) produces the solution to the FEM formulation of the scalar PDE problem

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

on a 2-D or 3-D region Ω, or the system PDE problem

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f,$

with geometry, mesh, and boundary conditions specified in model, with initial value u0 and initial derivative with respect to time ut0. The variables c, a, f, and d in the equation correspond to the function coefficients c, a, f, and d respectively.

example

u = hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d) solves the problem using boundary conditions b and finite element mesh specified in [p,e,t].

example

u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M) solves the problem based on finite element matrices that encode the equation, mesh, and boundary conditions.

u = hyperbolic(___,rtol) and u = hyperbolic(___,rtol,atol) modify the solution process by passing to the ODE solver a relative tolerance rtol, and optionally an absolute tolerance atol.

example

u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M,___,'DampingMatrix',D) modifies the problem to include a damping matrix D.

u = hyperbolic(___,'Stats','off') turns off the display of internal ODE solver statistics during the solution process.

## Examples

collapse all

### Hyperbolic Equation

Solve the wave equation

on the square domain specified by squareg.

Create a PDE model and import the geometry.

model = createpde;
geometryFromEdges(model,@squareg);
pdegplot(model,'EdgeLabels','on')
ylim([-1.1,1.1])
axis equal

Set Dirichlet boundary conditions for , and Neumann boundary conditions

for . (The Neumann boundary condition is the default condition, so the second specification is redundant.)

applyBoundaryCondition(model,'Edge',[2,4],'u',0);
applyBoundaryCondition(model,'Edge',[1,3],'g',0);

Set the initial conditions

u0 = 'atan(cos(pi/2*x))';
ut0 = '3*sin(pi*x).*exp(cos(pi*y))';

Set the solution times.

tlist = linspace(0,5,31);

Give coefficients for the problem.

c = 1;
a = 0;
f = 0;
d = 1;

Generate a mesh and solve the PDE.

generateMesh(model,'Hmax',0.1);
u1 = hyperbolic(u0,ut0,tlist,model,c,a,f,d);
549 successful steps
69 failed attempts
1238 function evaluations
1 partial derivatives
172 LU decompositions
1237 solutions of linear systems

Plot the solution at the first and last times.

figure
pdeplot(model,'xydata',u1(:,1))
figure
pdeplot(model,'xydata',u1(:,end))

For a version of this example with animation, run pdedemo6.

### Hyperbolic Equation using Legacy Syntax

Solve the wave equation

on the square domain specified by squareg, using a geometry function to specify the geometry, a boundary function to specify the boundary conditions, and using initmesh to create the finite element mesh.

Specify the geometry as @squareg and plot the geometry.

g = @squareg;
pdegplot(g,'EdgeLabels','on')
ylim([-1.1,1.1])
axis equal

Set Dirichlet boundary conditions for , and Neumann boundary conditions

for . (The Neumann boundary condition is the default condition, so the second specification is redundant.)

The squareb3 function specifies these boundary conditions.

b = @squareb3;

Set the initial conditions

u0 = 'atan(cos(pi/2*x))';
ut0 = '3*sin(pi*x).*exp(cos(pi*y))';

Set the solution times.

tlist = linspace(0,5,31);

Give coefficients for the problem.

c = 1;
a = 0;
f = 0;
d = 1;

Create a mesh and solve the PDE.

[p,e,t] = initmesh(g);
u = hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);
462 successful steps
70 failed attempts
1066 function evaluations
1 partial derivatives
156 LU decompositions
1065 solutions of linear systems

Plot the solution at the first and last times.

figure
pdeplot(p,e,t,'xydata',u(:,1))
figure
pdeplot(p,e,t,'xydata',u(:,end))

For a version of this example with animation, run pdedemo6.

### Hyperbolic Solution Using Finite Element Matrices

Solve a hyperbolic problem using finite element matrices.

Create a model and import the 'BracketWithHole.stl' geometry.

model = createpde();
importGeometry(model,'BracketWithHole.stl');
h = pdegplot(model,'FaceLabels','on');
h(1).FaceAlpha = 0.5;

Set coefficients c = 1, a = 0, f = 0.5, and d = 1.

c = 1;
a = 0;
f = 0.5;
d = 1;

Generate a mesh for the model.

generateMesh(model);

Create initial conditions and boundary conditions. The boundary condition for the rear face (face 3) is Dirichlet with value 0. All other faces have the default boundary condition. The initial condition is u(0) = 0, du/dt(0) = x/2. Give the initial condition on the derivative by calculating the x-position of each node in xpts, and passing x/2.

applyBoundaryCondition(model,'Face',3,'u',0);
u0 = 0;
xpts = model.Mesh.Nodes(1,:);
ut0 = xpts(:)/2;

Create the associated finite element matrices.

[Kc,Fc,B,ud] = assempde(model,c,a,f);
[~,M,~] = assema(model,0,d,f);

Solve the PDE for times from 0 to 2.

tlist = linspace(0,5,50);
u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M);
1486 successful steps
64 failed attempts
2981 function evaluations
1 partial derivatives
278 LU decompositions
2980 solutions of linear systems

View the solution at a few times. Scale all the plots to have the same color range by using the caxis command.

umax = max(max(u));
umin = min(min(u));

subplot(2,2,1)
pdeplot3D(model,'colormapdata',u(:,5));
caxis([umin umax]);
title('Time 1/2')
subplot(2,2,2)
pdeplot3D(model,'colormapdata',u(:,10));
caxis([umin umax]);
title('Time 1')
subplot(2,2,3)
pdeplot3D(model,'colormapdata',u(:,15));
caxis([umin umax]);
title('Time 3/2')
subplot(2,2,4)
pdeplot3D(model,'colormapdata',u(:,20));
caxis([umin umax]);
title('Time 2')

The solution seems to have a frequency of one, because the plots at times 1/2 and 3/2 show maximum values, and those at times 1 and 2 show minimum values.

### Hyperbolic Equation with Damping

Solve a hyperbolic problem that includes damping. You must use the finite element matrix form to use damping.

Create a model and import the 'BracketWithHole.stl' geometry.

model = createpde();
importGeometry(model,'BracketWithHole.stl');
h = pdegplot(model,'FaceLabels','on');
h(1).FaceAlpha = 0.5;

Set coefficients c = 1, a = 0, f = 0.5, and d = 1.

c = 1;
a = 0;
f = 0.5;
d = 1;

Generate a mesh for the model.

generateMesh(model);

Create initial conditions and boundary conditions. The boundary condition for the rear face (face 3) is Dirichlet with value 0. All other faces have the default boundary condition. The initial condition is u(0) = 0, du/dt(0) = x/2. Give the initial condition on the derivative by calculating the x-position of each node in xpts, and passing x/2.

applyBoundaryCondition(model,'Face',3,'u',0);
u0 = 0;
xpts = model.Mesh.Nodes(1,:);
ut0 = xpts(:)/2;

Create the associated finite element matrices.

[Kc,Fc,B,ud] = assempde(model,c,a,f);
[~,M,~] = assema(model,0,d,f);

Use a damping matrix that is 10% of the mass matrix.

Damping = 0.1*M;

Solve the PDE for times from 0 to 2.

tlist = linspace(0,5,50);
u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M,'DampingMatrix',Damping);
1419 successful steps
65 failed attempts
2745 function evaluations
1 partial derivatives
268 LU decompositions
2744 solutions of linear systems

Plot the maximum value at each time.

plot(max(u))
xlabel('Time')
ylabel('Maximum value')
title('Maximum of Solution')

The oscillations damp slightly as time increases.

## Input Arguments

collapse all

### u0 — Initial conditionvector | text expression

Initial condition, specified as a scalar, vector of nodal values, or text expression. The initial condition is the value of the solution u at the initial time, specified as a column vector of values at the nodes. The nodes are either p in the [p,e,t] data structure, or are model.Mesh.Nodes. For details, see Initial Conditions.

• If the initial condition is a constant scalar v, specify u0 as v.

• If there are Np nodes in the mesh, and N equations in the system of PDEs, specify u0 as a column vector of Np*N elements, where the first Np elements correspond to the first component of the solution u, the second Np elements correspond to the second component of the solution u, etc.

• Give a text expression of a function, such as 'x.^2 + 5*cos(x.*y)'. If you have a system of N > 1 equations, give a text array such as

char('x.^2 + 5*cos(x.*y)',...
'tanh(x.*y)./(1+z.^2)')

Example: x.^2+5*cos(y.*x)

Data Types: double | char
Complex Number Support: Yes

### ut0 — Initial derivativevector | text expression

Initial derivative, specified as a vector or text expression. The initial gradient is the value of the derivative of the solution u at the initial time, specified as a vector of values at the nodes. The nodes are either p in the [p,e,t] data structure, or are model.Mesh.Nodes. See Initial Conditions.

• If the initial derivative is a constant value v, specify u0 as v.

• If there are Np nodes in the mesh, and N equations in the system of PDEs, specify ut0 as a vector of Np*N elements, where the first Np elements correspond to the first component of the solution u, the second Np elements correspond to the second component of the solution u, etc.

• Give a text expression of a function, such as 'x.^2 + 5*cos(x.*y)'. If you have a system of N > 1 equations, use a text array such as

char('x.^2 + 5*cos(x.*y)',...
'tanh(x.*y)./(1+z.^2)')

For details, see Initial Conditions.

Example: p(1,:).^2+5*cos(p(2,:).*p(1,:))

Data Types: double | char
Complex Number Support: Yes

### tlist — Solution timesreal vector

Solution times, specified as a real vector. The solver returns the solution to the PDE at the solution times.

Example: 0:0.2:4

Data Types: double

### model — PDE modelPDEModel object

PDE model, specified as a PDEModel object.

Example: model = createpde(1)

### c — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. c represents the c coefficient in the scalar PDE

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specifyc in various ways, detailed in c Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: 'cosh(x+y.^2)'

Data Types: double | char | function_handle
Complex Number Support: Yes

### a — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. a represents the a coefficient in the scalar PDE

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

There are a wide variety of ways of specifying a, detailed in a or d Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: 2*eye(3)

Data Types: double | char | function_handle
Complex Number Support: Yes

### f — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. f represents the f coefficient in the scalar PDE

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specifyf in various ways, detailed in f Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: char('sin(x)';'cos(y)';'tan(z)')

Data Types: double | char | function_handle
Complex Number Support: Yes

### d — PDE coefficientscalar or matrix | character array | coefficient function

PDE coefficient, specified as a scalar or matrix, as a character array, or as a coefficient function. d represents the d coefficient in the scalar PDE

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\nabla u\right)+au=f,$

or in the system of PDEs

$d\frac{{\partial }^{2}u}{\partial {t}^{2}}-\nabla \cdot \left(c\otimes \nabla u\right)+au=f.$

You can specifyd in various ways, detailed in a or d Coefficient for Systems. See also Scalar PDE Coefficients, Specify Scalar PDE Coefficients in String Form, Specify 2-D Scalar Coefficients in Function Form, Specify 3-D PDE Coefficients in Function Form, and Coefficients for Systems of PDEs.

Example: 2*eye(3)

Data Types: double | char | function_handle
Complex Number Support: Yes

### b — Boundary conditionsboundary matrix | boundary file

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a string naming the file.

For more information on boundary conditions, see Forms of Boundary Condition Specification.

Example: b = 'circleb1' or equivalently b = @circleb1

Data Types: double | char | function_handle

### p — Mesh nodesoutput of initmesh | output of meshToPet

Mesh nodes, specified as the output of initmesh or meshToPet. For the structure of a p matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

### e — Mesh edgesoutput of initmesh | output of meshToPet

Mesh edges, specified as the output of initmesh or meshToPet. For the structure of e, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

### t — Mesh elementsoutput of initmesh | output of meshToPet

Mesh elements, specified as the output of initmesh or meshToPet. Mesh elements are the triangles or tetrahedra that form the finite element mesh. For the structure of a t matrix, see Mesh Data for [p,e,t] Triples: 2-D and Mesh Data for [p,e,t] Triples: 3-D.

Example: [p,e,t] = initmesh(g)

Data Types: double

### Kc — Stiffness matrixsparse matrix | full matrix

Stiffness matrix, specified as a sparse matrix or as a full matrix. See Elliptic Equations. Typically, Kc is the output of assempde.

Load vector, specified as a vector. See Elliptic Equations. Typically, Fc is the output of assempde.

### B — Dirichlet nullspacesparse matrix

Dirichlet nullspace, returned as a sparse matrix. See Algorithms. Typically, B is the output of assempde.

### ud — Dirichlet vectorvector

Dirichlet vector, returned as a vector. See Algorithms. Typically, ud is the output of assempde.

### M — Mass matrixsparse matrix | full matrix

Mass matrix. specified as a sparse matrix or a full matrix. See Elliptic Equations.

To obtain the input matrices for pdeeig, hyperbolic or parabolic, run both assema and assempde:

[Kc,Fc,B,ud] = assempde(model,c,a,f);
[~,M,~] = assema(model,0,d,f);
 Note:   Create the M matrix using assema with d, not a, as the argument before f.

Data Types: double
Complex Number Support: Yes

### rtol — Relative tolerance for ODE solver1e-3 (default) | positive real

Relative tolerance for ODE solver, specified as a positive real.

Example: 2e-4

Data Types: double

### atol — Absolute tolerance for ODE solver1e-6 (default) | positive real

Absolute tolerance for ODE solver, specified as a positive real.

Example: 2e-7

Data Types: double

### D — Damping matrixmatrix

Damping matrix, specified as a matrix. D has the same size as the stiffness matrix Kc or the mass matrix M. When you include D, hyperbolic solves the following ODE for the variable v:

${B}^{T}MB\frac{{d}^{2}v}{d{t}^{2}}+{B}^{T}DB\frac{dv}{dt}+Kv=F,$

with initial condition u0 and initial derivative ut0. Then hyperbolic returns the solution u = B*v + ud.

For an example using D, see Dynamics of a Damped Cantilever Beam.

Example: alpha*M + beta*K

Data Types: double
Complex Number Support: Yes

### Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example:

### 'Stats' — Display ODE solver statistics'on' (default) | 'off'

Display ODE solver statistics, specified as 'on' or 'off'. Suppress the display by setting Stats to 'off'.

Example: x = hyperbolic(u0,ut0,tlist,model,c,a,f,d,'Stats','off')

Data Types: char

## Output Arguments

collapse all

### u — PDE solutionmatrix

PDE solution, returned as a matrix. The matrix is Np*N-by-T, where Np is the number of nodes in the mesh, N is the number of equations in the PDE (N = 1 for a scalar PDE), and T is the number of solution times, meaning the number of elements of tlist. The solution matrix has the following structure.

• The first Np elements of each column in u represent the solution of equation 1, then next Np elements represent the solution of equation 2, etc. The solution u is the value at the corresponding node in the mesh.

• Column i of u represents the solution at time tlist(i).

To obtain the solution at an arbitrary point in the geometry, use pdeInterpolant.

To plot the solution, use pdeplot for 2-D geometry, or see Plot 3-D Solutions.

expand all

### Algorithms

hyperbolic internally calls assema, assemb, and assempde to create finite element matrices corresponding to the problem. It calls ode15s to solve the resulting system of ordinary differential equations. For details, see Hyperbolic Equations.