Ctrl+K

Дискретизация Задачи Коши. Частный случай, общий подход

Дискретизация Задачи Коши. Частный случай, общий подход#

Рассмотрим общий подход к численному решению ОДУ.

Задача Коши ставится следующим образом:

\[\begin{split} \left\{\begin{array}{l} \dfrac{d u}{d t}=f(t, u), \quad t \in [0, T] \\ u(0)=u_{0} \end{array}\right. \end{split}\]

Далее будем полагать, что решение существует и единственно.

В общем случае аналитически не решается - а решать-то хочется.

Понятное дело, на компьютере не можем работать с континуумом значений. Необходимо произвести дискретизацию задачи.

Зададим сетку с равным шагом:

\[ 0=t_0<t_1<t_2<\cdots<t_{n-1}<t_N=T \]
\[ \tau=t_{n+1}-t_n \]

Также зададим дискретную функцию \(y_n\) на этой сетке, от которой хотим

\[ y_n \approx u\left(t_n\right), \quad n=0,1, \cdots, N \]

Таким образом, \(y_n\) - это значения решения дискретизованного ОДУ на заданной сетке.

Важно. Заметьте, \(u(t)\) - истинное (математическое) решение ОДУ, а \(y_n\) - это сеточно-заданная функция. Таки разные вещи.

Вспомним про численную аппроксимацию производной. Выберем какую-нибудь численную схему, например, самую простую. Тогда

\[ \left.\frac{d u}{d t}\right|_{t_n} =\lim _{\tau \rightarrow 0} \frac{u(t_n+\tau)-u(t_n)}{\tau} \approx \frac{y_{n+1}-y_n}{\tau} \]

Подставляя это в изначальное уравнение, мы получаем дискретизованную Задачу Коши или разностная задача:

\[\begin{split} \left\{\begin{array}{l} \cfrac{y_{n+1}-y_n}{\tau}=f(t_n, y_n), \quad t \in [0, T] \\ y_0=u_{0} \end{array}\right. \end{split}\]

Примечание. Конкретно такой вид дискретизации называется явная схема Эйлера. Естественно, дискретизовывать можно по разному, выбирая различные разностные схемы и правые части, которые формально стремятся к истинной задаче Коши при \(\tau \rightarrow 0\).

Явную схему Эйлера уже понятно, как решать. Это же просто реккурента:

\[ y_{n+1}=y_n+\tau \cdot f\left(t_n, y_n\right), \quad n \geqslant 1 \]

Далее постараемся обобщить этот подход.

Hide code cell source 
from bokeh.io import output_notebook
from bokeh.resources import INLINE
output_notebook(INLINE)
Loading BokehJS ...
Hide code cell source 
import numpy as np
import sympy
from bokeh.io import show, output_file
from bokeh.layouts import column
from bokeh.models import ColumnDataSource, Slider, CustomJS, MultiLine
from bokeh.plotting import figure
from bokeh.palettes import Sunset10
blue, yellow = Sunset10[0], Sunset10[5]
from bokeh.models import Label

# Define the ODE to be solved
t = sympy.symbols('t')
y = sympy.Function('y')(t)
f = -y + sympy.sin(t)

# Set initial conditions and parameters
y0 = 1.0  # initial value of y
t0 = 0.0  # initial time
T = 10.0  # final time

# Calculate the exact solution using sympy and convert to numpy function
y_exact = sympy.dsolve(sympy.Eq(y.diff(t), f), y, ics={y.subs(t, t0): y0})
y_exact_func = sympy.lambdify(t, y_exact.rhs, 'numpy')
t_arr_dense = np.linspace(t0, T, 1000)  # denser time array for exact solution
y_exact_arr = y_exact_func(t_arr_dense)

# Define the initial value of h for the plot
h_init = 0.5

# Initialize arrays to store the results
t_arr = np.arange(t0, T + h_init, h_init)  # time array
y_arr = np.zeros_like(t_arr)  # solution array
y_arr[0] = y0  # set the initial value of y

# Solve the ODE using Explicit Euler scheme
for i in range(0, t_arr.size - 1):
    y_arr[i + 1] = y_arr[i] + h_init * f.subs([(y, y_arr[i]), (t, t_arr[i])]).evalf()

# Create a separate ColumnDataSource for the dense arrays
dense_source = ColumnDataSource(data=dict(t_arr_dense=t_arr_dense, y_exact_arr=y_exact_arr))

# Create the ColumnDataSource for the other arrays
source = ColumnDataSource(data=dict(
    t_arr=t_arr,
    y_arr=y_arr,
    y_exact_func=y_exact_func(t_arr),
))

# Create a separate ColumnDataSource for the vertical intervals
intervals_source = ColumnDataSource(data=dict(xs=[[t_arr[i], t_arr[i]] for i in range(len(t_arr))],
                                              ys=[[y_arr[i], y_exact_func(t_arr[i])] for i in range(len(t_arr))]))

# Create the plot
plot = figure(x_range=(t0, T), y_range=(-1.5, 1.1), width=500, height=300,
              title='ODE Solution with Explicit Euler Scheme')
plot.toolbar.active_drag = None
plot.line('t_arr_dense', 'y_exact_arr', line_color='black', source=dense_source)
plot.circle('t_arr', 'y_exact_func', radius = 0.05, line_color='red', fill_color=yellow, source=source)
plot.circle('t_arr', 'y_arr', radius = 0.05, line_color='blue', fill_color=blue, source=source)
plot.multi_line('xs', 'ys', line_color='grey', line_dash='dotted', source=intervals_source).level = 'underlay'  # Add vertical intervals

import sys
import os
try:
    sys.stderr = open(os.devnull, 'w') # for wsl build
except:
    pass

label = Label(
    text=r"$$y(x)\text{ - True Solution}$$",
    x=10, y=65,    
    x_units="screen", y_units="screen",
)


label1 = Label(
    text=r'$$y^{(τ)}=\{y_n\}_{n=0}^N \text{ - Discretized ODE Solution}$$',
    x=20, y=35,    
    x_units="screen", y_units="screen",
)
plot.circle(0.25, -1.05, radius = 0.05, line_color='blue', fill_color=blue)

label2 = Label(
    text=r"$$[y(x)]^{(τ)}=\{y(x_n)\}_{n=0}^N\text{ - True Solution Projection}$$",
    x=20, y=5,    
    x_units="screen", y_units="screen",
)
plot.circle(0.25, -1.35, radius = 0.05, line_color='red', fill_color=yellow)
plot.add_layout(label)
plot.add_layout(label1)
plot.add_layout(label2)

# Create a slider widget for h
slider = Slider(start=0.01, end=1.0, value=h_init, step=0.01, title='h', sizing_mode = "stretch_width")

# Define the JavaScript callback function
update_plot_code = """
const h = this.value;
const t0 = 0.0;
const T = 10.0;
const y0 = 1.0;
const t_arr = Array.from(Array(Math.ceil((T - t0) / h) + 1), (_, i) => t0 + i * h);
const y_arr = new Array(t_arr.length).fill(0);
y_arr[0] = y0;
for (let i = 0; i < t_arr.length - 1; i++) {
    const t = t_arr[i];
    const y = y_arr[i];
    y_arr[i + 1] = y + h * (Math.sin(t) - y);
}
const t_arr_dense = Array.from(Array(1000), (_, i) => t0 + i * (T - t0) / 999);
const y_exact_func = (t) => 1.5*Math.exp(-t) + Math.sin(t)/2 - Math.cos(t)/2;
const y_exact_arr = t_arr_dense.map(y_exact_func);
const y_exact_func_arr = t_arr.map(y_exact_func);

// Add code to create vertical intervals
const xs = [];
const ys = [];
for (let i = 0; i < t_arr.length; i++) {
    const y_min = Math.min(y_arr[i], y_exact_func_arr[i]);
    const y_max = Math.max(y_arr[i], y_exact_func_arr[i]);
    xs.push([t_arr[i], t_arr[i]]);
    ys.push([y_min, y_max]);
}

source.data = {
    t_arr: t_arr,
    y_arr: y_arr,
    y_exact_func: y_exact_func_arr,
};
dense_source.data = {
    t_arr_dense: t_arr_dense,
    y_exact_arr: y_exact_arr,
};
intervals_source.data = {
    xs: xs,
    ys: ys,
};
"""

callback = CustomJS(args=dict(source=source, dense_source=dense_source, intervals_source=intervals_source), code=update_plot_code)
slider.js_on_change('value', callback)

# Trigger the initial update of the plot
callback.code.replace("this.value", str(h_init))

# Create a layout with the plot and slider
layout = column(plot, slider)

# Display the layout
show(layout)

Примечание. Видим, что при больших значениях шага, решение дискретизованного ОДУ (синие точки) не совпадает с истинным (математическим) решением. Обратите так же внимание на обозначения в легенде.

Примеры различных простых дискретизаций#

Пусть задача Коши \(\dfrac{d y}{d x}=f(t, y)\)

Тогда, несколько вариантов соответствющих ей разностных задач \(\hat{L}_{\tau} y^{(\tau)}=F^{(\tau)}\):

  • Явная схема Эйлера

\[ \frac{y_{n+1}-y_n}{\tau}=f\left(t_n, \, y_n\right) \]

  • Неявная схема Эйлера

\[ \frac{y_{n+1}-y_n}{\tau}=f\left(t_{n+1}, \, y_{n+1}\right) \]

Неявная она потому, что в общем случае на каждом шаге приходится решать нелинейное уравнение на \(y_{n+1}\).

  • Симметризованная схема

\[ \frac{y_{n+1}-y_n}{\tau}=\frac{1}{2} \cdot \left(f\left(t_n, \, y_n\right)+f\left(t_{n+1}, \, y_{n+1}\right)\right) \]

  • Схема с центральной разностью

\[ \frac{x_{n+1}-x_{n-1}}{2 \cdot \tau}=f\left(x_{n}, t_{n}\right) \]

Видно, что строить их достаточно просто - главное обеспечить «аля» сходимость \(\hat{L}_{\tau} \rightarrow \hat{L}\) и \(F^{(\tau)} \rightarrow F(x)\).

Примечание. Первые три схемы относятся к двуслойным схемам, а четвёртая к трёхслойным. Подумайте почему.

Примечание. Начальные условия в виде производной в какой-то точке также необходимо дискретизовывать, и это влияет на точность. Смотреть пример ниже.

Для определения того, какая схема лучше и корректна ли схема вообще, необходимо ввести какие-то числовые параметры.

Содержание