计算科学 - 积分曲线的计算 - 吾爱随笔录

积分曲线的计算

计算科学颂

2021-11-30 03:58:47

考虑单位圆 $x^2 + y^2 = 1$ , 它也代表 ODE 的解

x + y \frac{d y}{d x} = 0; y (0) = 1

$x + y \frac{dy}{dx} = 0; y(0) = 1$

假设我们不知道如何解析求解上述方程（有更复杂的情况，解析解是不可能的），有限差分离散化给出：

y^{n + 1} = y^{n} - \frac{x^{n}}{y^{n}} δ, x^{0} = 0, y^{0} = 1.

$y^{n+1} = y^n - \frac{x^n}{y^n} \delta, x^0 = 0, y^0 = 1.$

δ

$\delta$ 是步长。该方法存在奇异性

(x = 1, y = 0)

$(x=1,y=0)$ ，其中圆的斜率发散到无穷大。

我想知道是否有更好的方法可以恢复整个圆圈。一个关键方面似乎是整合方向，即 $\delta$ , 必须在某个时候反转。

1个回答

我认为您可以将方程式重写为

x + \frac{1}{2} \frac{d u}{d x} = 0,

$x + \frac{1}{2}\frac{d u}{d x} = 0\, ,$

和 $u = y^2$ . 那么，差分方程将是

u^{n + 1} = u^{n} - 2 x Δ x .

$u^{n+1} = u^{n} - 2x \Delta x\, .$

可以通过计算结果的平方根来恢复原始变量。

通过这些更改，我设法在以下代码段中解决了它

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

#%% Solution
npts = 100
x = np.linspace(0, 1, npts)
u = np.ones_like(x)
dx = x[1] - x[0]
for cont in range(1, npts):
    u[cont] = u[cont - 1] - 2*x[cont - 1]*dx

#%% Plotting
plt.figure(figsize=(2, 2))
plt.plot(x, np.sqrt(u))
plt.axis("image")
plt.xlabel("x")
plt.ylabel("y")
plt.xticks([0, 0.5, 1])
plt.yticks([0, 0.5, 1])
plt.savefig("ode_sol.png")

而且，这就是情节

其它你可能感兴趣的问题

上一篇切比雪夫的网格映射下一篇最小化内存使用的 Eratosthenes 筛