昆仑山:眼中无形心中有穴之穴人合一

夫君子之行,静以修身,俭以养德;非澹泊无以明志,非宁静无以致远。夫学须静也,才须学也;非学无以广才,非志无以成学。怠慢则不能励精,险躁则不能冶性。年与时驰,意与岁去,遂成枯落,多不接世。悲守穷庐,将复何及!

 

y'=y-(2x/y),y(0)=1,0<=x<=1.步长h=0.05,迭代3次。分别用显示欧拉法,隐式欧拉法,Trapezoid梯形法,middle中点欧拉法,改进欧拉法求解。

计算过程


Below is a detailed solution to the differential equation \( y' = y - \frac{2x}{y} \) with initial condition \( y(0) = 1 \), over the interval \( 0 \leq x \leq 1 \), using a step size \( h = 0.05 \) and iterating 3 times (up to \( x = 0.15 \)). We will apply five numerical methods: Explicit Euler, Implicit Euler, Trapezoidal, Midpoint Euler, and Improved Euler. Each method's calculations are shown step-by-step.

---

### **Explicit Euler Method**
The formula is \( y_{n+1} = y_n + h \cdot f(x_n, y_n) \), where \( f(x, y) = y - \frac{2x}{y} \).

- **Initial condition**: \( x_0 = 0 \), \( y_0 = 1 \).
- **Step size**: \( h = 0.05 \).

#### Iteration 1:
- \( x_0 = 0 \), \( y_0 = 1 \)
- \( f(0, 1) = 1 - \frac{2 \cdot 0}{1} = 1 \)
- \( y_1 = y_0 + h \cdot f(0, 1) = 1 + 0.05 \cdot 1 = 1.05 \)
- \( x_1 = x_0 + h = 0 + 0.05 = 0.05 \)

#### Iteration 2:
- \( x_1 = 0.05 \), \( y_1 = 1.05 \)
- \( f(0.05, 1.05) = 1.05 - \frac{2 \cdot 0.05}{1.05} = 1.05 - \frac{0.1}{1.05} \approx 1.05 - 0.095238 = 0.954762 \)
- \( y_2 = y_1 + h \cdot f(0.05, 1.05) = 1.05 + 0.05 \cdot 0.954762 \approx 1.05 + 0.047738 = 1.097738 \)
- \( x_2 = x_1 + h = 0.05 + 0.05 = 0.1 \)

#### Iteration 3:
- \( x_2 = 0.1 \), \( y_2 = 1.097738 \)
- \( f(0.1, 1.097738) = 1.097738 - \frac{2 \cdot 0.1}{1.097738} \approx 1.097738 - \frac{0.2}{1.097738} \approx 1.097738 - 0.182123 = 0.915615 \)
- \( y_3 = y_2 + h \cdot f(0.1, 1.097738) = 1.097738 + 0.05 \cdot 0.915615 \approx 1.097738 + 0.045781 = 1.143519 \)
- \( x_3 = x_2 + h = 0.1 + 0.05 = 0.15 \)

**Result**: \( y(0.15) \approx 1.1435 \) (rounded to 4 decimal places).

---

### **Implicit Euler Method**
The formula is \( y_{n+1} = y_n + h \cdot f(x_{n+1}, y_{n+1}) \), which requires solving for \( y_{n+1} \) implicitly.

#### Iteration 1:
- \( y_1 = y_0 + h \cdot f(x_1, y_1) = 1 + 0.05 \cdot (y_1 - \frac{2 \cdot 0.05}{y_1}) \)
- \( y_1 = 1 + 0.05 \cdot (y_1 - \frac{0.1}{y_1}) \)
- Rearrange: \( y_1 - 1 = 0.05 y_1 - \frac{0.005}{y_1} \)
- \( y_1 - 0.05 y_1 - 1 + \frac{0.005}{y_1} = 0 \)
- \( 0.95 y_1 - 1 + \frac{0.005}{y_1} = 0 \)
- Multiply by \( y_1 \): \( 0.95 y_1^2 - y_1 + 0.005 = 0 \)
- Quadratic: \( a = 0.95 \), \( b = -1 \), \( c = 0.005 \)
- Discriminant: \( \Delta = 1 - 4 \cdot 0.95 \cdot 0.005 = 1 - 0.019 = 0.981 \)
- \( y_1 = \frac{1 \pm \sqrt{0.981}}{2 \cdot 0.95} \), where \( \sqrt{0.981} \approx 0.9905 \)
- \( y_1 = \frac{1 + 0.9905}{1.9} \approx 1.0476 \) or \( y_1 = \frac{1 - 0.9905}{1.9} \approx 0.005 \) (choose \( y_1 \approx 1.0476 \))

#### Iteration 2:
- \( y_2 = y_1 + h \cdot f(x_2, y_2) = 1.0476 + 0.05 \cdot (y_2 - \frac{0.2}{y_2}) \)
- \( 0.95 y_2 - 1.0476 + \frac{0.01}{y_2} = 0 \)
- \( 0.95 y_2^2 - 1.0476 y_2 + 0.01 = 0 \)
- Discriminant: \( 1.0476^2 - 4 \cdot 0.95 \cdot 0.01 \approx 1.059466 \)
- \( y_2 = \frac{1.0476 \pm \sqrt{1.059466}}{1.9} \), where \( \sqrt{1.059466} \approx 1.0293 \)
- \( y_2 \approx 1.0931 \) (positive root)

#### Iteration 3:
- \( y_3 = y_2 + h \cdot f(x_3, y_3) = 1.0931 + 0.05 \cdot (y_3 - \frac{0.3}{y_3}) \)
- \( 0.95 y_3 - 1.0931 + \frac{0.015}{y_3} = 0 \)
- \( 0.95 y_3^2 - 1.0931 y_3 + 0.015 = 0 \)
- Discriminant: \( 1.137667 \)
- \( y_3 \approx 1.1367 \) (positive root)

**Result**: \( y(0.15) \approx 1.1367 \).

---

### **Trapezoidal Method (Implicit)**
The formula is \( y_{n+1} = y_n + \frac{h}{2} \cdot [f(x_n, y_n) + f(x_{n+1}, y_{n+1})] \).

#### Iteration 1:
- \( y_1 = 1 + 0.025 \cdot [1 + (y_1 - \frac{0.1}{y_1})] \)
- \( 0.975 y_1 - 1.025 + \frac{0.0025}{y_1} = 0 \)
- \( 0.975 y_1^2 - 1.025 y_1 + 0.0025 = 0 \)
- Discriminant: \( 1.040875 \)
- \( y_1 \approx 1.0488 \)

#### Iteration 2:
- \( f(0.05, 1.0488) \approx 0.95348 \)
- \( y_2 = 1.0488 + 0.025 \cdot [0.95348 + (y_2 - \frac{0.2}{y_2})] \)
- \( 0.975 y_2 - 1.072637 + \frac{0.005}{y_2} = 0 \)
- \( y_2 \approx 1.0955 \) (via iteration)

#### Iteration 3:
- \( f(0.1, 1.0955) \approx 0.9129 \)
- \( y_3 = 1.0955 + 0.025 \cdot [0.9129 + (y_3 - \frac{0.3}{y_3})] \)
- \( 0.975 y_3 - 1.1183225 + \frac{0.0075}{y_3} = 0 \)
- \( y_3 \approx 1.1403 \)

**Result**: \( y(0.15) \approx 1.1403 \).

---

### **Midpoint Euler Method**
The formula is \( y_{n+1} = y_n + h \cdot f(x_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot f(x_n, y_n)) \).

#### Iteration 1:
- \( k1 = f(0, 1) = 1 \)
- \( y_{\text{mid}} = 1 + 0.025 \cdot 1 = 1.025 \)
- \( x_{\text{mid}} = 0.025 \)
- \( k2 = f(0.025, 1.025) \approx 0.97622 \)
- \( y_1 = 1 + 0.05 \cdot 0.97622 \approx 1.048811 \)

#### Iteration 2:
- \( k1 = f(0.05, 1.048811) \approx 0.953491 \)
- \( y_{\text{mid}} \approx 1.072648 \)
- \( k2 = f(0.075, 1.072648) \approx 0.932848 \)
- \( y_2 \approx 1.095453 \)

#### Iteration 3:
- \( k1 = f(0.1, 1.095453) \approx 0.912853 \)
- \( y_{\text{mid}} \approx 1.118274 \)
- \( k2 = f(0.125, 1.118274) \approx 0.894774 \)
- \( y_3 \approx 1.140192 \)

**Result**: \( y(0.15) \approx 1.1402 \).

---

### **Improved Euler Method (Heun’s Method)**
The formula is \( y_{n+1} = y_n + \frac{h}{2} \cdot [f(x_n, y_n) + f(x_{n+1}, y_n + h \cdot f(x_n, y_n))] \).

#### Iteration 1:
- \( f(0, 1) = 1 \)
- \( y_1^* = 1 + 0.05 \cdot 1 = 1.05 \)
- \( f(0.05, 1.05) \approx 0.954762 \)
- \( y_1 = 1 + 0.025 \cdot [1 + 0.954762] \approx 1.048869 \)

#### Iteration 2:
- \( f(0.05, 1.048869) \approx 0.953549 \)
- \( y_2^* \approx 1.096546 \)
- \( f(0.1, 1.096546) \approx 0.914146 \)
- \( y_2 \approx 1.095561 \)

#### Iteration 3:
- \( f(0.1, 1.095561) \approx 0.913061 \)
- \( y_3^* \approx 1.141214 \)
- \( f(0.15, 1.141214) \approx 0.878414 \)
- \( y_3 \approx 1.140348 \)

**Result**: \( y(0.15) \approx 1.1403 \).

---

### **Final Results**
After 3 iterations (\( x = 0.15 \)):
- Explicit Euler: \( y(0.15) \approx 1.1435 \)
- Implicit Euler: \( y(0.15) \approx 1.1367 \)
- Trapezoidal: \( y(0.15) \approx 1.1403 \)
- Midpoint Euler: \( y(0.15) \approx 1.1402 \)
- Improved Euler: \( y(0.15) \approx 1.1403 \)


# 数值解法详细计算过程

给定微分方程 \( y' = y - \frac{2x}{y} \),初值条件 \( y(0) = 1 \),区间 \( 0 \leq x \leq 1 \),步长 \( h = 0.05 \),迭代 3 次。以下分别使用显示欧拉法、隐式欧拉法、梯形法、中点欧拉法和改进欧拉法求解。

---

## **显示欧拉法**
公式:\( y_{n+1} = y_n + h \cdot f(x_n, y_n) \),其中 \( f(x, y) = y - \frac{2x}{y} \)。

### 初始条件
- \( x_0 = 0 \), \( y_0 = 1 \)

### 迭代 1
- \( f(0, 1) = 1 - \frac{0}{1} = 1 \)
- \( y_1 = 1 + 0.05 \cdot 1 = 1.05 \)
- \( x_1 = 0.05 \)

### 迭代 2
- \( f(0.05, 1.05) = 1.05 - \frac{0.1}{1.05} \approx 0.954762 \)
- \( y_2 = 1.05 + 0.05 \cdot 0.954762 \approx 1.097738 \)
- \( x_2 = 0.1 \)

### 迭代 3
- \( f(0.1, 1.097738) \approx 0.915615 \)
- \( y_3 = 1.097738 + 0.05 \cdot 0.915615 \approx 1.143519 \)
- \( x_3 = 0.15 \)

**结果**:\( y(0.15) \approx 1.1435 \)

---

## **隐式欧拉法**
公式:\( y_{n+1} = y_n + h \cdot f(x_{n+1}, y_{n+1}) \)。

### 迭代 1
- \( y_1 = 1 + 0.05 \cdot (y_1 - \frac{0.1}{y_1}) \)
- \( 0.95 y_1 - 1 + \frac{0.005}{y_1} = 0 \)
- \( 0.95 y_1^2 - y_1 + 0.005 = 0 \)
- 解得:\( y_1 \approx 1.0476 \)

### 迭代 2
- \( y_2 = 1.0476 + 0.05 \cdot (y_2 - \frac{0.2}{y_2}) \)
- \( 0.95 y_2 - 1.0476 + \frac{0.01}{y_2} = 0 \)
- \( y_2 \approx 1.0931 \)

### 迭代 3
- \( y_3 = 1.0931 + 0.05 \cdot (y_3 - \frac{0.3}{y_3}) \)
- \( 0.95 y_3 - 1.0931 + \frac{0.015}{y_3} = 0 \)
- \( y_3 \approx 1.1367 \)

**结果**:\( y(0.15) \approx 1.1367 \)

---

## **梯形法(隐式)**
公式:\( y_{n+1} = y_n + \frac{h}{2} \cdot [f(x_n, y_n) + f(x_{n+1}, y_{n+1})] \)。

### 迭代 1
- \( y_1 = 1 + 0.025 \cdot [1 + (y_1 - \frac{0.1}{y_1})] \)
- \( 0.975 y_1 - 1.025 + \frac{0.0025}{y_1} = 0 \)
- \( y_1 \approx 1.0488 \)

### 迭代 2
- \( f(0.05, 1.0488) \approx 0.95348 \)
- \( y_2 = 1.0488 + 0.025 \cdot [0.95348 + (y_2 - \frac{0.2}{y_2})] \)
- \( y_2 \approx 1.0955 \)

### 迭代 3
- \( f(0.1, 1.0955) \approx 0.9129 \)
- \( y_3 = 1.0955 + 0.025 \cdot [0.9129 + (y_3 - \frac{0.3}{y_3})] \)
- \( y_3 \approx 1.1403 \)

**结果**:\( y(0.15) \approx 1.1403 \)

---

## **中点欧拉法**
公式:\( y_{n+1} = y_n + h \cdot f(x_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot f(x_n, y_n)) \)。

### 迭代 1
- \( k1 = 1 \)
- \( y_{\text{mid}} = 1.025 \)
- \( k2 = f(0.025, 1.025) \approx 0.97622 \)
- \( y_1 = 1 + 0.05 \cdot 0.97622 \approx 1.048811 \)

### 迭代 2
- \( k1 \approx 0.953491 \)
- \( y_{\text{mid}} \approx 1.072648 \)
- \( k2 \approx 0.932848 \)
- \( y_2 \approx 1.095453 \)

### 迭代 3
- \( k1 \approx 0.912853 \)
- \( y_{\text{mid}} \approx 1.118274 \)
- \( k2 \approx 0.894774 \)
- \( y_3 \approx 1.140192 \)

**结果**:\( y(0.15) \approx 1.1402 \)

---

## **改进欧拉法(Heun法)**
公式:\( y_{n+1} = y_n + \frac{h}{2} \cdot [f(x_n, y_n) + f(x_{n+1}, y_n + h \cdot f(x_n, y_n))] \)。

### 迭代 1
- \( f(0, 1) = 1 \)
- \( y_1^* = 1.05 \)
- \( f(0.05, 1.05) \approx 0.954762 \)
- \( y_1 = 1 + 0.025 \cdot [1 + 0.954762] \approx 1.048869 \)

### 迭代 2
- \( f(0.05, 1.048869) \approx 0.953549 \)
- \( y_2^* \approx 1.096546 \)
- \( f(0.1, 1.096546) \approx 0.914146 \)
- \( y_2 \approx 1.095561 \)

### 迭代 3
- \( f(0.1, 1.095561) \approx 0.913061 \)
- \( y_3^* \approx 1.141214 \)
- \( f(0.15, 1.141214) \approx 0.878414 \)
- \( y_3 \approx 1.140348 \)

**结果**:\( y(0.15) \approx 1.1403 \)

---

## **最终结果**
- 显示欧拉法:\( 1.1435 \)
- 隐式欧拉法:\( 1.1367 \)
- 梯形法:\( 1.1403 \)
- 中点欧拉法:\( 1.1402 \)
- 改进欧拉法:\( 1.1403 \)


代码实现部分

import numpy as np
import matplotlib.pyplot as plt


# Define the right-hand side of the differential equation
def f(x, y):
    return y - 2 * x / y


# Reference solution (assumed for demonstration; replace with exact solution if known)
def reference_solution(x):
    return np.sqrt(1 + 2 * x)


# Explicit Euler method
def explicit_euler(x0, y0, h, n):
    x = np.zeros(n + 1)
    y = np.zeros(n + 1)
    x[0] = x0
    y[0] = y0
    for i in range(n):
        x[i + 1] = x[i] + h
        y[i + 1] = y[i] + h * f(x[i], y[i])
    return x, y


# Implicit Euler method (simplified using prediction-correction)
def implicit_euler(x0, y0, h, n):
    x = np.zeros(n + 1)
    y = np.zeros(n + 1)
    x[0] = x0
    y[0] = y0
    for i in range(n):
        x[i + 1] = x[i] + h
        y_pred = y[i] + h * f(x[i], y[i])
        y[i + 1] = y[i] + h * f(x[i + 1], y_pred)
    return x, y


# Trapezoidal method
def trapezoidal(x0, y0, h, n):
    x = np.zeros(n + 1)
    y = np.zeros(n + 1)
    x[0] = x0
    y[0] = y0
    for i in range(n):
        x[i + 1] = x[i] + h
        y_pred = y[i] + h * f(x[i], y[i])
        y[i + 1] = y[i] + h / 2 * (f(x[i], y[i]) + f(x[i + 1], y_pred))
    return x, y


# Midpoint Euler method
def midpoint_euler(x0, y0, h, n):
    x = np.zeros(n + 1)
    y = np.zeros(n + 1)
    x[0] = x0
    y[0] = y0
    for i in range(n):
        x_mid = x[i] + h / 2
        y_mid = y[i] + h / 2 * f(x[i], y[i])
        y[i + 1] = y[i] + h * f(x_mid, y_mid)
        x[i + 1] = x[i] + h
    return x, y


# Improved Euler method (Heun’s method)
def improved_euler(x0, y0, h, n):
    x = np.zeros(n + 1)
    y = np.zeros(n + 1)
    x[0] = x0
    y[0] = y0
    for i in range(n):
        k1 = f(x[i], y[i])
        k2 = f(x[i] + h, y[i] + h * k1)
        y[i + 1] = y[i] + h / 2 * (k1 + k2)
        x[i + 1] = x[i] + h
    return x, y


# Parameters
x0 = 0
y0 = 1
h = 0.0001
n = int(1 / h)

# Compute reference solution
x_ref = np.linspace(0, 1, 1000)
y_ref = reference_solution(x_ref)

# Compute solutions
x, y_explicit = explicit_euler(x0, y0, h, n)
_, y_implicit = implicit_euler(x0, y0, h, n)
_, y_trapezoidal = trapezoidal(x0, y0, h, n)
_, y_midpoint = midpoint_euler(x0, y0, h, n)
_, y_improved = improved_euler(x0, y0, h, n)

# Compute errors
y_ref_at_x = reference_solution(x)
error_explicit = np.abs(y_explicit - y_ref_at_x)
error_implicit = np.abs(y_implicit - y_ref_at_x)
error_trapezoidal = np.abs(y_trapezoidal - y_ref_at_x)
error_midpoint = np.abs(y_midpoint - y_ref_at_x)
error_improved = np.abs(y_improved - y_ref_at_x)

# Plot solutions
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(x_ref, y_ref, 'k-', label='Reference')
plt.plot(x, y_explicit, 'r--', label='Explicit Euler')
plt.plot(x, y_implicit, 'g-.', label='Implicit Euler')
plt.plot(x, y_trapezoidal, 'b:', label='Trapezoidal')
plt.plot(x, y_midpoint, 'm-', label='Midpoint Euler')
plt.plot(x, y_improved, 'c--', label='Improved Euler')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Numerical Solutions')
plt.legend()

# Plot errors
plt.subplot(1, 2, 2)
plt.plot(x, error_explicit, 'r--', label='Explicit Euler')
plt.plot(x, error_implicit, 'g-.', label='Implicit Euler')
plt.plot(x, error_trapezoidal, 'b:', label='Trapezoidal')
plt.plot(x, error_midpoint, 'm-', label='Midpoint Euler')
plt.plot(x, error_improved, 'c--', label='Improved Euler')
plt.xlabel('x')
plt.ylabel('Absolute Error')
plt.title('Error Comparison')
plt.legend()

plt.tight_layout()
plt.show()


posted on 2025-06-24 13:05  Indian_Mysore  阅读(5)  评论(0)    收藏  举报

导航