Wolfram Language & System Documentation Center

InitializePDECoefficients

NDSolve`FEM`

InitializePDECoefficients

InitializePDECoefficients[vd,sd,rules]

initializes the coefficients specified by rules in accordance with variable data vd and solution data sd to generate a PDECoefficientData object.

Details and Options

The coefficients are assumed to come from a second-order system of PDEs in space dimensions of the form:

In InitializePDECoefficients[vd,sd,rules], the rules should be of the form name->coefficient. The possible coefficient names are:

"LoadCoefficients"	{{f₁},{f₂},…}	is a scalar
"LoadDerivativeCoefficients"	{{γ₁},{γ₂},…}	is a vector of length
"DiffusionCoefficients"	{{-c₁₁,-c₁₂,…},{-c₂₁,-c₂₂,…},…}	may be specified as scalars, diagonal vectors of length , or × matrices
"ConservativeConvectionCoefficients"	{{-α₁₁,-α₁₂,…},{-α₂₁,-α₂₂,…},…}	is a vector of length
"ConvectionCoefficients"	{{β₁₁,β₁₂,…},{β₂₁,β₂₂,…},…}	is a vector of length
"ReactionCoefficients"	{{a₁₁,a₁₂,…},{a₂₁,a₂₂,…},…}	is a scalar

If a rule is not specified for any of these coefficient names, the coefficients of that type are all assumed to be 0.
DiscretizePDE will take the generated PDECoefficientData object and provide system matrices such that

where is the stiffness matrix and a load vector.
Transient systems of PDEs may be specified up through second, with the first-order temporal derivative and the second-order temporal derivative.
The transient coefficient form is:

DiscretizePDE will take the generated PDECoefficientData object and provide system matrices such that

where is the mass matrix, the damping matrix, the stiffness matrix and a load vector.
The possible damping coefficient names are:

"DampingDiffusionCoefficients"	{{-dc₁₁,-dc₁₂,…},{-dc₂₁,-dc₂₂,…},…}	may be specified as scalars, diagonal vectors of length or × matrices
"DampingConservativeConvectionCoefficients"	{{-dα₁₁,-dα₁₂,…},{-dα₂₁,-dα₂₂,…},…}	is a vector of length
"DampingConvectionCoefficients"	{{dβ₁₁,dβ₁₂,…},{dβ₂₁,dβ₂₂,…},…}	is a vector of length
"DampingReactionCoefficients"	{{dα₁₁,dα₁₂,…},{dα₂₁,dα₂₂,…},…}	is a scalar

The possible mass coefficient names are:

"MassDiffusionCoefficients"	{{-mc₁₁,-mc₁₂,…},{-mc₂₁,-mc₂₂,…},…}	may be specified as scalars, diagonal vectors of length or × matrices
"MassConservativeConvectionCoefficients"	{{-mα₁₁,-mα₁₂,…},{-mα₂₁,-mα₂₂,…},…}	is a vector of length
"MassConvectionCoefficients"	{{mβ₁₁,mβ₁₂,…},{mβ₂₁,mβ₂₂,…},…}	is a vector of length
"MassReactionCoefficients"	{{mα₁₁,mα₁₂,…},{mα₂₁,mα₂₂,…},…}	is a scalar

"Mass" coefficients contribute to the mass matrix , "Damping" coefficients contribute to the damping matrix , and the remaining coefficients contribute to the stiffness matrix and the load vector .
NDSolve reduces transient systems so that they are first order in time.
The coefficients can be functions of space, time, parameters, dependent variables and first-order derivatives of dependent variables.
Variable data vd and solution data sd are corresponding lists of variables and values. Templates for vd and sd may be generated using NDSolve`VariableData and NDSolve`SolutionData, and components may be set using NDSolve`SetSolutionDataComponent.
InitializePDECoefficients verifies and optimizes the coefficients in accordance with variable data vd and solution data sd.
The "Space" component of vd and sd should be set to the spatial variables and the spatial mesh represented as a NumericalRegion object, respectively.
The "DependentVariables" component of vd should be a list of dependent variable name symbols without arguments.
For time-dependent problems, the "Time" component of vd and sd should be set to the temporal variable and the initial time, respectively.
For nonlinear problems, the "DependentVariables" component of sd should be set to the initial seedings for the dependent variables.
For parametric problems, the "Parameters" component of vd and sd should be set to the parametric variables and the initial parametric values, respectively.
InitializePDECoefficients has the following options:
"VerificationData" Automatic specify PDE coefficient verification data
Options given to InitializePDECoefficients can be given to NDSolve by specifying "InitializePDECoefficientsOptions". »
Setting the option from NDSolve and related functions is explained in NDSolve Finite Element Options.

Examples

open all close all

Basic Examples (2)

Load the finite element package:

Wolfram Language code: Needs["NDSolve`FEM`"]

Set up a NumericalRegion:

Wolfram Language code: nr = ToNumericalRegion[Disk[]]

Set up variable and solution data:

Wolfram Language code:

vd = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> nr];

Convert a Laplace equation into coefficients:

Wolfram Language code: InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}]

Convert a Poisson equation into coefficients:

Wolfram Language code: InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, "LoadCoefficients" -> {{1}}]

Scope (6)

Coefficients can be functions:

Wolfram Language code: f[x_, y_] := {{x ^ 2, 0}, {x + y, y ^ 2}}

Wolfram Language code: InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-f[x, y]}}]

Coefficients can be interpolation functions:

Wolfram Language code: ifun = Interpolation[Flatten[Table[Evaluate[{{x, y}, Sin[x y]}], {x, 0, 1, 0.2}, {y, 0, 1 / 2, 0.1}], 1]]

Wolfram Language code: InitializePDECoefficients[vd, sd, "LoadCoefficients" -> {{-ifun[x, y]}}]

Coefficients can be nonlinear:

Wolfram Language code:

sdnonlinear = sd;
NDSolve`SetSolutionDataComponent[sdnonlinear, "DependentVariables", {0.}];
InitializePDECoefficients[vd, sdnonlinear, "DiffusionCoefficients" -> {{{{-u[x, y], 0}, {0, -Derivative[1, 0][u][x, y]}}}}]

If the input for "DiffusionCoefficient" is a vector, it is interpreted as the diagonal of the requested input matrix:

Wolfram Language code:

InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-{1, 1}}}] === InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}]

If the input for "DiffusionCoefficient" is a scalar, the scalar is put in all places of the diagonal of the requested matrix:

Wolfram Language code:

InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-1}}] === InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}]

Initialize a system of PDEs:

Wolfram Language code:

vd2 = vd;
NDSolve`SetSolutionDataComponent[vd2, "DependentVariables", {u, v}];
InitializePDECoefficients[vd2, sd, "DiffusionCoefficients" -> -{{{{(100000/91), 0}, {0, (5000/13)}}, {{0, (30000/91)}, {(5000/13), 0}}}, {{{0, (5000/13)}, {(30000/91), 0}}, {{(5000/13), 0}, {0, (100000/91)}}}}]

Options (1)

"VerificationData" (1)

The PDE coefficients that are given to InitializePDECoefficients are verified so that they evaluate to numeric input:

Wolfram Language code: rc[x_ ? NumericQ, y_ ? NumericQ] := (Print["Test Coordinate: ", {x, y}];1)

Wolfram Language code: InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{-rc[x, y]}}]

The verification is done by finding a test coordinate that lies within the NumericalRegion. It may, however, happen that the automatically chosen test coordinate or test ElementMarker are at a singularity:

Wolfram Language code: rc[x_ ? NumericQ, y_ ? NumericQ] := (Print["Test Coordinate: ", {x, y}];1 / (x + 0.5 + y + 0.5))

Wolfram Language code: InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{-rc[x, y]}}]

Overwriting the automatic finding of the verification data can be done with the "VerificationData" option:

Wolfram Language code:

InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{-rc[x, y]}}, "VerificationData" -> {"Coordinate" -> {0, 0}}]

Additionally, the automatic ElementMarker can be changed:

Wolfram Language code:

InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{-rc[x, y]}}, "VerificationData" -> {"Coordinate" -> {0, 0}, "ElementMarker" -> 1}]

For nonlinear problems, the verification data for the dependents and the derivative of the dependents can be modified with "DependentVariables" and "DependentVariablesDerivative", respectively:

Wolfram Language code:

sdnonlinear = sd;
NDSolve`SetSolutionDataComponent[sdnonlinear, "DependentVariables", {0.}];
InitializePDECoefficients[vd, sdnonlinear, "ReactionCoefficients" -> {{-1 / u[x, y]}}]

Change the test value for the dependent variable to be 1 and not 0 from the initial condition:

Wolfram Language code:

InitializePDECoefficients[vd, sdnonlinear, "ReactionCoefficients" -> {{-1 / u[x, y]}}, "VerificationData" -> {"DependentVariables" -> {1}}]

Change the test value of the derivative of the dependent variable to the vector {1,1}. The components of the vector are the values of the derivative in each space direction. Currently, a default 0.1 is used:

Wolfram Language code:

InitializePDECoefficients[vd, sdnonlinear, "ReactionCoefficients" -> {{-1 / (0.1 - Derivative[1, 0][u][x, y])}}, "VerificationData" -> {"DependentVariablesDerivatives" -> {{1, 1}}}]

Applications (6)

Specify a Laplace equation over a disk with Dirichlet and Neumann conditions for and for :

Wolfram Language code:

pded = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}];
bcd = InitializeBoundaryConditions[vd, sd, {{DirichletCondition[u[x, y] == 0, x ≤ -0.3], NeumannValue[1., x ≥ 0.35]}}];
md = InitializePDEMethodData[vd, sd];

Inspect the diffusion coefficients:

Wolfram Language code: pded["DiffusionCoefficients"]

Discretize, solve, interpolate and visualize the solution:

Wolfram Language code:

dpde = DiscretizePDE[pded, md, sd];
dbc = DiscretizeBoundaryConditions[bcd, md, sd];
{l, s, d, m} = dpde["SystemMatrices"];
DeployBoundaryConditions[{l, s, d, m}, dbc];
res = LinearSolve[s, l];
efun = ElementMeshInterpolation[{nr["ElementMesh"]}, res];
Plot3D[efun[x, y], {x, y}∈nr["ElementMesh"]]

The "LoadCoefficients" coefficient is used by NIntegrate to integrate over a region:

Wolfram Language code:

ipde = InitializePDECoefficients[vd, sd, "LoadCoefficients" -> {{1}}];
md = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[ipde, md, sd];
Total[dpde["LoadVector"], 2]

Use NIntegrate to verify the result:

Wolfram Language code: NIntegrate[1, {x, y}∈Disk[]]

A "ReactionCoefficients" of can be used to evaluate a function, for example, over a region:

Wolfram Language code:

ipde = InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{1}}, "LoadCoefficients" -> {{x ^ 2 + y ^ 2}}];
md = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[ipde, md, sd];
sol = LinearSolve[dpde["StiffnessMatrix"], dpde["LoadVector"]];
if = ElementMeshInterpolation[{md["ElementMesh"]}, sol]

Visualize the evaluated function:

Wolfram Language code: Plot3D[if[x, y], {x, y}∈Disk[]]

A "ReactionCoefficients" of in combination with a "LoadCoefficients" can be used to scale a function, for example, over a region by a factor of :

Wolfram Language code:

ipde = InitializePDECoefficients[vd, sd, "ReactionCoefficients" -> {{10}}, "LoadCoefficients" -> {{x ^ 2 + y ^ 2}}];
md = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[ipde, md, sd];
sol = LinearSolve[dpde["StiffnessMatrix"], dpde["LoadVector"]];
if = ElementMeshInterpolation[{md["ElementMesh"]}, sol]

Visualize the scaled evaluated function:

Wolfram Language code: Plot3D[if[x, y], {x, y}∈Disk[], PlotRange -> All]

A "ReactionCoefficients" in combination with a "LoadDerivativeCoefficients" can be used to interpolate the derivative of a function over a region:

Wolfram Language code:

uExact[x_] := Cos[x] + 1
vd1D = NDSolve`VariableData[{"DependentVariables", "Space"} -> {{u}, {x}}];
sd1D = NDSolve`SolutionData["Space" -> ToNumericalRegion[Line[{{0}, {π}}]]];
ipde = InitializePDECoefficients[vd1D, sd1D, "ReactionCoefficients" -> {{-1}}, 
	"LoadDerivativeCoefficients" -> {{{uExact[x]}}}];
ibcd = InitializeBoundaryConditions[vd1D, sd1D, {{DirichletCondition[u[x] == (D[uExact[x], x] /. x -> 0), x == 0], DirichletCondition[u[x] == (D[uExact[x], x] /. x -> π), x == π]}}];
md = InitializePDEMethodData[vd1D, sd1D];
dpde = DiscretizePDE[ipde, md, sd1D];
dbcd = DiscretizeBoundaryConditions[ibcd, md, sd1D, "Stationary"];
{s, l} = {dpde["StiffnessMatrix"], dpde["LoadVector"]};
DeployBoundaryConditions[{l, s}, dbcd];
sol = LinearSolve[s, l];
if = ElementMeshInterpolation[{md["ElementMesh"]}, sol]

Visualize the difference between the numerical approximation and the exact derivative:

Wolfram Language code: Plot[Evaluate[{if[x] - Div[{uExact[x]}, {x}]}], {x, 0, π}, PlotRange -> All]

Specify an equation with mixed spatial and temporal derivatives in the domain :

Wolfram Language code: aeqn = -Derivative[0, 2][V][t, x] - Derivative[1, 2][V][t, x]

Specify Dirichlet conditions , and initial conditions :

Wolfram Language code:

ic = V[0, x] == 0;
bcRight = V[t, 1] == 0;
bcLeftValue = Sin[2 * π * t];
bcLeft = V[t, 0] == bcLeftValue;
tEnd = 3;

Set up the variable and solution data:

Wolfram Language code:

vd1D = NDSolve`VariableData[{"DependentVariables" -> {u}, "Space" -> {x}, "Time" -> t}];
nr1D = ToNumericalRegion[FullRegion[1], {{0, 1}}];
sd1D = NDSolve`SolutionData[{"Space" -> nr1D, "Time" -> 0.}];

Initialize the coefficients:

Wolfram Language code:

initCoeffs = InitializePDECoefficients[vd1D, sd1D, 
	"DiffusionCoefficients" -> {{-IdentityMatrix[1]}}, 
	"DampingDiffusionCoefficients" -> {{-IdentityMatrix[1]}}
	]

Initialize the boundary conditions:

Wolfram Language code:

initBCs = InitializeBoundaryConditions[vd1D, sd1D, {{DirichletCondition[u[t, x] == bcLeftValue, x == 0], DirichletCondition[u[t, x] == 0, x == 1]}}]

Initialize the method data:

Wolfram Language code: md = InitializePDEMethodData[vd1D, sd1D]

Discretize the PDE coefficients. In this case, all PDE coefficients are time independent:

Wolfram Language code: sdpde = DiscretizePDE[initCoeffs, md, sd1D, "Stationary"];

Inspect the stiffness matrix:

Wolfram Language code: sdpde["StiffnessMatrix"]

Inspect the damping matrix:

Wolfram Language code: sdpde["DampingMatrix"]

Discretize the boundary conditions:

Wolfram Language code: sbcs = DiscretizeBoundaryConditions[initBCs, md, sd1D]

Write a function to compute the residual of the PDE:

Wolfram Language code:

discretizePDEResidual[t_ ? NumericQ, u_ ? VectorQ, dudt_ ? VectorQ] := 
	Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs}, 
	NDSolve`SetSolutionDataComponent[sd1D, "Time", t];
	NDSolve`SetSolutionDataComponent[sd1D, "DependentVariables", u];
	l = sdpde["LoadVector"];
	s = sdpde["StiffnessMatrix"];
	d = sdpde["DampingMatrix"];
	tbcs = DiscretizeBoundaryConditions[initBCs, md, sd1D, "Transient"];
	DeployBoundaryConditions[{l, s, d}, tbcs];
	DeployBoundaryConditions[{l, s, d}, sbcs];
	d.dudt + s.u - l
	]

Evaluate the initial conditions:

Wolfram Language code: init = EvaluateOnElementMesh[{x}, 0, md["ElementMesh"]]["ValuesOnGrid"];

Set up a function for the sparsity patterns of the system matrices:

Wolfram Language code:

sparsityPatterns[t_, u_ ? VectorQ] := Module[
	{s, d, tdpde, nldpde}, 
	NDSolve`SetSolutionDataComponent[sd1D, "Time", t];
	NDSolve`SetSolutionDataComponent[sd1D, "DependentVariables", u];
	s = sdpde["StiffnessMatrix"];
	d = sdpde["DampingMatrix"];
	{s["PatternArray"], d["PatternArray"]}
	]

Compute the sparsity patterns:

Wolfram Language code: {stiffnessSparsity, dampingSparsity} = sparsityPatterns[10 ^ -6, init]

Solve the equation:

Wolfram Language code:

ufun = NDSolveValue[{discretizePDEResidual[t, U[t], U'[t]] == 0, U[0] == init}, U, {t, 0, tEnd}, Method -> {"EquationSimplification" -> "Residual"}, Jacobian -> {"X" -> {Automatic, Sparse -> stiffnessSparsity}, "XP" -> {Automatic, Sparse -> dampingSparsity}}]

Create a unified interpolating function:

Wolfram Language code:

femSolution = ElementMeshInterpolation[{ufun["Coordinates"][[1]], md["ElementMesh"]}, Partition[ufun["ValuesOnGrid"], 1]]

For comparison, compute a solution with the tensor product grid method:

Wolfram Language code: tpgSolution = NDSolveValue[{aeqn == 0, bcLeft, bcRight, ic}, V, {x, 0, 1}, {t, 0, tEnd}]

Plot the difference between the TPG and FEM methods:

Wolfram Language code: Plot3D[tpgSolution[t, x] - femSolution[t, x], {t, 0, tEnd}, {x, 0., 1.}]

Properties & Relations (2)

Options given to InitializePDECoefficients can be given to NDSolve by specifying "InitializePDECoefficientOptions":

Wolfram Language code:

NDSolve[{Laplacian[u[x, y], {x, y}] == 1 / (x + 1 / 2 + y + 1 / 2), DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[]]

Specify "VerificationData" to pass a coordinate that is not at a singularity:

Wolfram Language code:

NDSolve[{Laplacian[u[x, y], {x, y}] == 1 / (x + 1 / 2 + y + 1 / 2), DirichletCondition[u[x, y] == 0, True]}, u, {x, y}∈Disk[], Method -> {"FiniteElement", "InitializePDECoefficientsOptions" -> {"VerificationData" -> {"Coordinate" -> {0, 0}}}}]

Convert a Poisson equation into coefficients:

Wolfram Language code: pdec = InitializePDECoefficients[vd, sd, "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, "LoadCoefficients" -> {{1}}]

Use GetInactivePDE to create an inactive version of the input:

Wolfram Language code: GetInactivePDE[pdec, vd]

Possible Issues (1)

Typically, NDSolve will detect if coefficients do not evaluate to numeric quantities with the correct dimensions and warn about that:

Wolfram Language code:

NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] == k Sin[x], u[0, x] == 0, DirichletCondition[u[t, x] == 0, True]}, u, {t, 0, π}, {x, -π, π}, Method -> {"PDEDiscretization" -> {"MethodOfLines"}}]

To verify that the coefficients given are adequate for each operator, the coefficients undergo a test evaluation using the initial point in time and a coordinate from the region. If the combination of the coefficient and the used values cancel out, NDSolve cannot detect if a coefficient has not been specified. NDSolve will then give a message that it was not able to compute the elements contributed by a specific operator:

Wolfram Language code:

NDSolveValue[{D[u[t, x], {t, 1}] - D[u[t, x], {x, 2}] == Cos[t] Sin[x] + k Sin[t] Sin[x], u[0, x] == 0, DirichletCondition[u[t, x] == 0, True]}, u, {t, 0, π}, {x, -π, π}, Method -> {"PDEDiscretization" -> {"MethodOfLines"}}]

In the preceding case, evaluating the load operator at t=0 and x=-π will let the load operator vanish, and NDSolve cannot detect that k has not been specified:

Wolfram Language code: (Cos[t] * Sin[x] + k * Sin[t] * Sin[x]) /. {t -> 0., x -> -π}

Top

More Learning

Tech Support

Wolfram Solutions

Wolfram Solutions For Education

Get Started

Grow Your Skills

Work with Us

Educational Programs for Adults

Educational Programs for Youth

Read

InitializePDECoefficients

Details and Options

Examples

Basic Examples (2)

Scope (6)

Options (1)

"VerificationData" (1)

Applications (6)

Properties & Relations (2)

Possible Issues (1)

Text

CMS

APA

BibTeX

BibLaTeX

InitializePDECoefficients

Details and Options

Examples

Basic Examples (2)

Scope (6)

Options (1)

"VerificationData" (1)

Applications (6)

Properties & Relations (2)

Possible Issues (1)

See Also

Tech Notes

Related Guides

Text

CMS

APA

BibTeX

BibLaTeX