# 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

```
[1]:
```

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

```
[2]:
```

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

Check `oph`

, *outer oriented port-Hamiltonian*, with

```
[3]:
```

```
oph.pr()
```

```
[3]:
```

```
<Figure size 800x600 with 1 Axes>
```

We can take out the knowns and label them by

```
[4]:
```

```
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\),

```
[5]:
```

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

.

```
[6]:
```

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

We now apply a particular discretization to this weak formulation,

```
[7]:
```

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

```
[8]:
```

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

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

```
[9]:
```

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

```
[9]:
```

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