```python
import sympy as sp
import pandas as pd
```

# Unequally spaced nodes


```python
# Helper function to compute divided differences
def div_diffs(nodes, values, max_order):
    if len(nodes) != len(values):
        raise ValueError("nodes and values lists must have the same length")
    diffs = [values]
    for order in range(max_order):
        current_diffs = []
        for i in range(1, len(nodes) - order):
            current_diffs.append((diffs[-1][i] - diffs[-1][i-1]) / (nodes[i + order] - nodes[i-1]))
        diffs.append(current_diffs)
    return diffs
```


```python
# Helper function to nicely print divided differences
def print_diffs(diffs):
    rows = 2 * len(diffs[0]) - 1
    for r in range(rows):
        row = '\t' if r % 2 == 1 else ''
        for i in range(r % 2, min(rows - r, min(r + 1, len(diffs))), 2):
            row += f'{float(diffs[i][(r - i) // 2]):.4f}\t\t'
        print(row[:-2])
```


```python
# Create DataFrame from the provided data with addition of comments
def create_usn_dataframe(new_node, closest_nodes, estimations, error, derivative_est, error_est):
    index = ['Closest nodes', f'P_i({new_node})', 'Actual error', 'M_(i+1)', f'R_i({new_node})']
    rows = [
        closest_nodes,
        estimations,
        error,
        derivative_est,
        error_est
    ]
    return pd.DataFrame(rows, index=index)
```


```python
x = sp.symbols('x')
f = sp.exp(x)
```


```python
# Nodes, function and new node are all taken from the book
nodes = [-0.3, -0.2, -0.1, 0, 0.1, 0.3]
new_node = 0.2
order = 3
values = [f.subs(x, n).evalf() for n in nodes]
```


```python
dds = div_diffs(nodes, values, order)
print_diffs(dds)
```

    0.7408
    	0.7791
    0.8187		0.4097
    	0.8611		0.1436
    0.9048		0.4528
    	0.9516		0.1587
    1.0000		0.5004
    	1.0517		0.1800
    1.1052		0.5724
    	1.2234
    1.3499



```python
for pos in range(len(nodes)):
    if nodes[pos] > new_node:
        break
# pos is where we'd put the new node among the others
```


```python
# Find position of the closest by value nodes
i, j = pos - 1, pos
pos_around_new = []
while (j - i) <= order + 1:
    if i < 0 and j >= len(nodes):
        break
    if j >= len(nodes) or new_node - nodes[i] <= nodes[j] - new_node:
        pos_around_new.append(i)
        i -= 1
    else:
        pos_around_new.append(j)
        j += 1
```


```python
nodes_around_new = [nodes[i] for i in pos_around_new]
```


```python
polynom = 0
mult = 1
estimated_values = []
for o in range(order + 1):
    m = min(pos_around_new[:o+1])
    polynom += dds[o][m] * mult
    estimated_values.append(polynom.subs(x, new_node).evalf())
    mult *= x - nodes[pos_around_new[o]]
print(f'Estimated value: {estimated_values[-1]}\nTrue value: {f.subs(x, new_node).evalf()}')
```

    Estimated value: 1.22143043351646
    True value: 1.22140275816017



```python
error = [f.subs(x, new_node) - p for p in estimated_values]
```


```python
# all the derivaties of e^x are e^x, and maximum of e^x on any segment [a, b], where it exists, is e^b.
# also, e^x is always > 0.
ms = []
for i in range(len(nodes_around_new)):
    ms.append(f.evalf(subs={x: max([new_node] + nodes_around_new[:i+1])}))
```


```python
# Compute estimated error
rs = []
mult = 1
for i in range(len(nodes_around_new)):
    mult *= abs(new_node - nodes_around_new[i]) / (i + 1)
    rs.append(ms[i] * mult)
```


```python
create_usn_dataframe(new_node, nodes_around_new, estimated_values, error, ms, rs)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Closest nodes</th>
      <td>0.3</td>
      <td>0.1</td>
      <td>0</td>
      <td>-0.1</td>
    </tr>
    <tr>
      <th>P_i(0.2)</th>
      <td>1.34985880757600</td>
      <td>1.22751486282583</td>
      <td>1.22179052060098</td>
      <td>1.22143043351646</td>
    </tr>
    <tr>
      <th>Actual error</th>
      <td>-0.128456049415833</td>
      <td>-0.00611210466565559</td>
      <td>-0.000387762440812178</td>
      <td>-2.76753562922227e-5</td>
    </tr>
    <tr>
      <th>M_(i+1)</th>
      <td>1.34985880757600</td>
      <td>1.34985880757600</td>
      <td>1.34985880757600</td>
      <td>1.34985880757600</td>
    </tr>
    <tr>
      <th>R_i(0.2)</th>
      <td>0.134985880757600</td>
      <td>0.00674929403788001</td>
      <td>0.000449952935858668</td>
      <td>3.37464701894001e-5</td>
    </tr>
  </tbody>
</table>
</div>



# Equally spaced nodes


```python
def finite_differences(nodes, values, max_order):
    if len(nodes) != len(values):
        raise ValueError("nodes and values lists must have the same length")
    diffs = [values]
    for order in range(max_order):
        current_diffs = []
        for i in range(1, len(nodes) - order):
            current_diffs.append(diffs[-1][i] - diffs[-1][i-1])
        diffs.append(current_diffs)
    return diffs
```


```python
def create_esn_dataframe(new_node, fin_diffs, ns, estimations, errors, err_est):
    index = ['∆k*y_(n-k)', 'N_k(t)', 'N_k*∆k*y_n(n-k)',
            f'P_k({new_node})', f'f({new_node}) - P_k({new_node})',
            '|R_k({new_node})| <=']
    rows = [fin_diffs, ns, [f * n for f, n in zip(fin_diffs, ns)],
           estimations, errors, err_est]
    return pd.DataFrame(rows, index=index)
```


```python
f = sp.cos(x)
a, b = -3.5, -2.5
new_node = -2.55
step = 0.1
order = 4
```


```python
nodes, values = [], []
c = a
while c <= b:
    nodes.append(round(c, 5))
    values.append(f.evalf(subs={x: c}))
    c = round(c + step, 5)
```


```python
for n, v in zip(nodes, values):
    print(f'f({n:.3g})\t= {v:.5g}')
```

    f(-3.5)	= -0.93646
    f(-3.4)	= -0.96680
    f(-3.3)	= -0.98748
    f(-3.2)	= -0.99829
    f(-3.1)	= -0.99914
    f(-3)	= -0.98999
    f(-2.9)	= -0.97096
    f(-2.8)	= -0.94222
    f(-2.7)	= -0.90407
    f(-2.6)	= -0.85689
    f(-2.5)	= -0.80114



```python
fds = finite_differences(nodes, values, order)
```


```python
print_diffs(fds)
```

    -0.9365
    	-0.0303
    -0.9668		0.0097
    	-0.0207		0.0002
    -0.9875		0.0099		-0.0001
    	-0.0108		0.0001
    -0.9983		0.0100		-0.0001
    	-0.0008		0.0000
    -0.9991		0.0100		-0.0001
    	0.0091		-0.0001
    -0.9900		0.0099		-0.0001
    	0.0190		-0.0002
    -0.9710		0.0097		-0.0001
    	0.0287		-0.0003
    -0.9422		0.0094		-0.0001
    	0.0382		-0.0004
    -0.9041		0.0090		-0.0001
    	0.0472		-0.0005
    -0.8569		0.0086
    	0.0557
    -0.8011



```python
# new_node comes at the end of the table, and we take that into account to find the polynom
b - step / 2 <= new_node and new_node < b
```




    True




```python
t = (x - b) / step
```


```python
# Compute N_k for k up `order` to estimate error for polynom of order `order - 1`
ns = [1]
for i in range(order + 1):
    ns.append(ns[-1] * (t + i) / (i + 1))
```


```python
# Calculate exact values of N_k in new_node
n_values = [ns[0]] + [n.evalf(subs={x: new_node}) for n in ns[1:]]
```


```python
# Calculate all the polynoms with order from 0 to `order`
polynoms = []
for i in range(order + 1):
    prev = polynoms[-1] if len(polynoms) > 0 else 0
    polynoms.append(prev + ns[i] * fds[i][-1])
```


```python
true_value = f.evalf(subs={x: new_node})
estimations = [p.evalf(subs={x: new_node}) for p in polynoms]
errors = [true_value - e for e in estimations]
```


```python
errors
```




    [-0.0289099196882884,
     -0.00103735077728173,
     3.28678694556661e-5,
     3.40281804211973e-6,
     -1.22832288784736e-7]




```python
turning_point = -sp.pi
df = sp.diff(f)
ms = []
for i in range(order + 1):
    p = b - (i+1) * step
    m = max(abs(df.evalf(subs={x: b})),
            abs(df.evalf(subs={x: p})))
    if p > turning_point:
        m = max(m, abs(df.evalf(subs={x: turning_point})))
    ms.append(m)
    
    df = sp.diff(df)
```


```python
# Error estimation
rs = [(ms[i] * abs(ns[i+1]) * step ** (i + 1)).evalf(subs={x: new_node}) for i in range(order + 1)]
```


```python
# Put all the finite differences values we used into a list
used_fds = [fds[i][-1] for i in range(order+1)]
```


```python
create_esn_dataframe(new_node, used_fds, n_values[:-1], estimations, errors, rs)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
      <th>3</th>
      <th>4</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>∆k*y_(n-k)</th>
      <td>-0.801143615546934</td>
      <td>0.0557451378220136</td>
      <td>0.00856174917389974</td>
      <td>-0.000471440822617297</td>
      <td>-9.02566484699552e-5</td>
    </tr>
    <tr>
      <th>N_k(t)</th>
      <td>1</td>
      <td>-0.499999999999998</td>
      <td>-0.125000000000000</td>
      <td>-0.0625000000000001</td>
      <td>-0.0390625000000001</td>
    </tr>
    <tr>
      <th>N_k*∆k*y_n(n-k)</th>
      <td>-0.801143615546934</td>
      <td>-0.0278725689110067</td>
      <td>-0.00107021864673747</td>
      <td>2.94650514135811e-5</td>
      <td>3.52565033085763e-6</td>
    </tr>
    <tr>
      <th>P_k(-2.55)</th>
      <td>-0.801143615546934</td>
      <td>-0.829016184457940</td>
      <td>-0.830086403104678</td>
      <td>-0.830056938053264</td>
      <td>-0.830053412402933</td>
    </tr>
    <tr>
      <th>f(-2.55) - P_k(-2.55)</th>
      <td>-0.0289099196882884</td>
      <td>-0.00103735077728173</td>
      <td>3.28678694556661e-5</td>
      <td>3.40281804211973e-6</td>
      <td>-1.22832288784736e-7</td>
    </tr>
    <tr>
      <th>|R_k({new_node})| &lt;=</th>
      <td>0.0299236072051977</td>
      <td>0.00125000000000000</td>
      <td>3.74045090064973e-5</td>
      <td>3.90625000000001e-6</td>
      <td>1.63644726903426e-7</td>
    </tr>
  </tbody>
</table>
</div>



f7c354ed7b14e53fa12b43dc7b9806faa7a8e93a
ipynb
Jupyter Notebook
year-2/computational-workshop/task-1.ipynb
Sergobot/university
7cd8c07fc660f1e19127c6488991ddd59d99643c
MIT
2019-09-05T08:43:47.000Z
2019-09-05T08:43:52.000Z
year-2/computational-workshop/task-1.ipynb
Sergobot/university
7cd8c07fc660f1e19127c6488991ddd59d99643c
MIT
year-2/computational-workshop/task-1.ipynb
Sergobot/university
7cd8c07fc660f1e19127c6488991ddd59d99643c
MIT
2020-10-04T07:40:20.000Z
2021-09-18T07:12:08.000Z
Qwen/Qwen-72B
1. YES
2. YES