Documentation

This is machine translation

Translated by
Mouseover text to see original. Click the button below to return to the English verison of the page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

hyperbolic

(Not recommended) 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.

`hyperbolic` is not recommended. Use `solvepde` instead.

Syntax

``u = hyperbolic(u0,ut0,tlist,model,c,a,f,d)``
``u = hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)``
``u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M)``
``u = hyperbolic(___,rtol)``
``````u = hyperbolic(___,rtol,atol)``````
``u = hyperbolic(u0,ut0,tlist,Kc,Fc,B,ud,M,___,'DampingMatrix',D)``
``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

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,'dirichlet','Edge',[2,4],'u',0); applyBoundaryCondition(model,'neumann','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,'GeometricOrder','linear','Hmax',0.1); u1 = hyperbolic(u0,ut0,tlist,model,c,a,f,d);```
```462 successful steps 51 failed attempts 1028 function evaluations 1 partial derivatives 135 LU decompositions 1027 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, see Wave Equation on a Square Domain.

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, see Wave Equation on a Square Domain.

Solve a hyperbolic problem using finite element matrices.

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

```model = createpde(); importGeometry(model,'BracketWithHole.stl'); figure pdegplot(model,'FaceLabels','on') view(30,30) title('Bracket with Face Labels')```

```figure pdegplot(model,'FaceLabels','on') view(-134,-32) title('Bracket with Face Labels, Rear View')```

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 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',4,'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);```
```1493 successful steps 70 failed attempts 2937 function evaluations 1 partial derivatives 276 LU decompositions 2936 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.

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'); figure pdegplot(model,'FaceLabels','on') view(30,30) title('Bracket with Face Labels')```

```figure pdegplot(model,'FaceLabels','on') view(-134,-32) title('Bracket with Face Labels, Rear View')```

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 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',4,'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);```
```1441 successful steps 70 failed attempts 2850 function evaluations 1 partial derivatives 288 LU decompositions 2849 solutions of linear systems ```

Plot the maximum value at each time. The oscillations damp slightly as time increases.

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

Input Arguments

collapse all

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 Solve PDEs with 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

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 Solve PDEs with 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 Solve PDEs with Initial Conditions.

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

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

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`

PDE model, specified as a `PDEModel` object.

Example: `model = createpde`

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 specify`c` in various ways, detailed in c Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

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

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

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 Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

Example: `2*eye(3)`

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

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 specify`f` in various ways, detailed in f Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

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

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

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 specify`d` in various ways, detailed in a or d Coefficient for Systems. See also Specify Scalar PDE Coefficients in Character Form, Specify 2-D Scalar Coefficients in Function Form, and Specify 3-D PDE Coefficients in Function Form.

Example: `2*eye(3)`

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

Boundary conditions, specified as a boundary matrix or boundary file. Pass a boundary file as a function handle or as a file name.

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

Data Types: `double` | `char` | `function_handle`

Mesh points, specified as a 2-by-`Np` matrix of points, where `Np` is the number of points in the mesh. For a description of the (`p`,`e`,`t`) matrices, see Mesh Data.

Typically, you use the `p`, `e`, and `t` data exported from the PDE Modeler app, or generated by `initmesh` or `refinemesh`.

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

Data Types: `double`

Mesh edges, specified as a `7`-by-`Ne` matrix of edges, where `Ne` is the number of edges in the mesh. For a description of the (`p`,`e`,`t`) matrices, see Mesh Data.

Typically, you use the `p`, `e`, and `t` data exported from the PDE Modeler app, or generated by `initmesh` or `refinemesh`.

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

Data Types: `double`

Mesh triangles, specified as a `4`-by-`Nt` matrix of triangles, where `Nt` is the number of triangles in the mesh. For a description of the (`p`,`e`,`t`) matrices, see Mesh Data.

Typically, you use the `p`, `e`, and `t` data exported from the PDE Modeler app, or generated by `initmesh` or `refinemesh`.

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

Data Types: `double`

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`.

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

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

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

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

Example: `2e-4`

Data Types: `double`

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

Example: `2e-7`

Data Types: `double`

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','off'`

collapse all

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

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 length 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 and Their Gradients.

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.