# How to discretize equations#

Now we demonstrate how to discretize a PDE object.

Pre-coded sample objects are stored in sample attribute of phyem. Invoke these samples by

:

import sys
ph_dir = '../../../../../'   # the path to dir that containing the phyem package.
sys.path.append(ph_dir)
import phyem as ph  # import the phyem package
ph.config._set_matplot_block(False)
samples = ph.samples


The partial differential equations of the canocical linear por-Hamiltonian are pre-coded. Call it through

:

oph = samples.pde_canonical_pH(n=3, p=3, periodic=True)  # where o on oph means outer


Check oph, outer oriented port-Hamiltonian, with

:

oph.pr()

:

<Figure size 800x600 with 1 Axes>


We can take out the knowns and label them by

:

a3, b2 = oph.unknowns


We now test oph with test functions from the spaces where a3 and b2 come from, and label the test functions by $$v^3$$ and $$u^2$$,

:

wf = oph.test_with(oph.unknowns, sym_repr=[r'v^3', r'u^2'])


Now, we apply integration by parts to the term indexed by '1-1'.

:

wf = wf.derive.integration_by_parts('1-1')


We now apply a particular discretization to this weak formulation,

:

td = wf.td
td.set_time_sequence()  # initialize a time sequence
td.define_abstract_time_instants('k-1', 'k-1/2', 'k')
td.differentiate('0-0', 'k-1', 'k')
td.average('0-1', b2, ['k-1', 'k'])

td.differentiate('1-0', 'k-1', 'k')
td.average('1-1', a3, ['k-1', 'k'])
dt = td.time_sequence.make_time_interval('k-1', 'k')

wf = td()

wf.unknowns = [
a3 @ td.time_sequence['k'],
b2 @ td.time_sequence['k'],
]

wf = wf.derive.split(
'0-0', 'f0',
[a3 @ td.ts['k'], a3 @ td.ts['k-1']],
['+', '-'],
factors=[1/dt, 1/dt],
)

wf = wf.derive.split(
'0-2', 'f0',
[ph.d(b2 @ td.ts['k-1']), ph.d(b2 @ td.ts['k'])],
['+', '+'],
factors=[1/2, 1/2],
)

wf = wf.derive.split(
'1-0', 'f0',
[b2 @ td.ts['k'], b2 @ td.ts['k-1']],
['+', '-'],
factors=[1/dt, 1/dt]
)

wf = wf.derive.split(
'1-2', 'f0',
[a3 @ td.ts['k-1'], a3 @ td.ts['k']],
['+', '+'],
factors=[1/2, 1/2],
)

wf = wf.derive.rearrange(
{
0: '0, 3 = 2, 1',
1: '3, 0 = 2, 1',
}
)


We now can write the weak formulation with matrix proxies. Before doing that, we need to

:

ph.space.finite(3)
mp = wf.mp()


The matrix format of the weak formulation leads to a linear system,

:

ls = mp.ls()
ls.pr()

:

<Figure size 1200x600 with 1 Axes>


Note that, till this moment, everything is still abstract. To do the numerical simulation, we need to bring them to a particular implementation, for example, the msepy, standing for mimetic spectral elements python, by calling ph.fem.apply function, which will be shown in other notebooks.