你的这段 Fortran 程序实现了 **拉格朗日多项式插值（Lagrange Interpolation）**，用于对函数 $ f(x) = \frac{1}{1 + x^2} $ 在区间 $[-5, 5]$ 上进行数值逼近，并分析误差。
这个函数是著名的 **“龙格函数”（Runge's function）**，常被用来展示高次多项式插值的 **龙格现象（Runge Phenomenon）** —— 即在端点附近出现剧烈振荡、误差增大。
---
## ✅ 代码结构总览
```fortran
program lagrange_reconstruction
implicit none
! 变量声明与精度定义
! 输入节点数 N
! 分配并生成等距插值节点 x_nodes 和函数值 y_nodes
! 在更密集的点上计算精确值和插值结果，输出对比表
! 包含：原函数 f(x)，拉格朗日插值子程序
end program
```
下面我们逐部分详细解释。
---
### 🔹 1. 精度设置：双精度浮点
```fortran
integer, parameter :: dp = selected_real_kind(15, 300)
```
- `dp` 表示 "double precision"
- `selected_real_kind(15, 300)` 要求：
- 至少 15 位有效数字（精度）
- 指数范围至少 ±300
- 这确保了使用机器支持的最高精度（通常是 64 位 double）
> ✅ 推荐做法：避免使用 `real*8`，这是非标准语法；用 `real(dp)` 更可移植。
---
### 🔹 2. 主变量声明
```fortran
integer :: N, i, j, num_points
real(dp), dimension(:), allocatable :: x_nodes, y_nodes
real(dp) :: a, b, h, x
real(dp) :: exact, interpolated, error
```
- `N`: 插值节点个数（从 0 到 N，共 N+1 个点）
- `x_nodes`, `y_nodes`: 存储插值节点及其函数值
- `a`, `b`: 区间边界 [-5, 5]
- `h`: 输出采样步长
- `x`: 当前评估点
- `exact`: 真实函数值 $ f(x) $
- `interpolated`: 拉格朗日插值近似值
- `error`: 绝对误差 $ |f(x) - p_N(x)| $
---
### 🔹 3. 设置区间和读入 N
```fortran
a = -5.0_dp
b = 5.0_dp
write(*,*) '请输入插值节点数 N :'
read(*,*) N
```
注意：这里 `N` 是最大索引，所以实际有 `N+1` 个节点！
例如你输入 `N=10`，就有 11 个等距节点。
---
### 🔹 4. 分配内存并生成数据
```fortran
allocate(x_nodes(0:N), y_nodes(0:N))
do i = 0, N
x_nodes(i) = a + (b - a) * i / N
y_nodes(i) = f(x_nodes(i))
end do
```
- 使用等距节点：$ x_i = -5 + \frac{10i}{N},\quad i=0,\dots,N $
- 计算每个节点上的函数值 $ y_i = f(x_i) = \frac{1}{1+x_i^2} $
⚠️ 注意：这些是 **插值所用的数据点**，之后将构造一个次数为 N 的多项式通过所有这些点。
---
### 🔹 5. 输出标题
```fortran
write(*,'(A)') 'x, exact_f(x), interpolated, error'
```
打印表头，后续每行输出四列数据。
---
### 🔹 6. 在密集网格上评估插值效果
```fortran
num_points = 101
h = (b - a) / (num_points - 1)
do i = 0, num_points - 1
x = a + i * h
exact = f(x)
call lagrange_interpolation(N, x_nodes, y_nodes, x, interpolated)
error = abs(exact - interpolated)
write(*,'(F6.2, 3E14.6)') x, exact, interpolated, error
end do
```
- 在 `[a,b]` 上取 101 个均匀点来检验插值质量
- 对每个点 `x`：
- 计算真实值 `exact = f(x)`
- 调用 `lagrange_interpolation` 得到插值近似值
- 计算绝对误差
- 输出格式：
- `F6.2`: 小数点前最多 4 位，保留 2 位小数（如 `-5.00`）
- `3E14.6`: 三个科学计数法数，宽度 14，6 位小数
---
### 🔹 7. 函数 `f(x)` 定义（Runge 函数）
```fortran
function f(x) result(res)
real(dp), intent(in) :: x
real(dp) :: res
res = 1.0_dp / (1.0_dp + x*x)
end function f
```
- 标准 Runge 函数：光滑但解析延拓在复平面有极点（±i），导致实轴上高阶插值不稳定。
- 高次插值时会出现严重 **龙格现象**（两端剧烈震荡）
---
### 🔹 8. 子程序：拉格朗日插值核心算法
```fortran
subroutine lagrange_interpolation(N, x_nodes, y_nodes, x, interpolated_value)
integer, intent(in) :: N
real(dp), dimension(0:N), intent(in) :: x_nodes, y_nodes
real(dp), intent(in) :: x
real(dp), intent(out) :: interpolated_value
integer :: i, j
real(dp) :: product, term
interpolated_value = 0.0_dp
do i = 0, N
product = 1.0_dp
do j = 0, N
if (j /= i) then
product = product * (x - x_nodes(j)) / (x_nodes(i) - x_nodes(j))
end if
end do
term = y_nodes(i) * product
interpolated_value = interpolated_value + term
end do
end subroutine
```
#### 📌 原理回顾：拉格朗日插值公式
$$
p_N(x) = \sum_{i=0}^{N} y_i \cdot \ell_i(x),\quad
\ell_i(x) = \prod_{\substack{j=0\\j\ne i}}^{N} \frac{x - x_j}{x_i - x_j}
$$
- 外层循环 `i`：遍历每个基函数 $\ell_i(x)$
- 内层循环 `j`：计算乘积项
- `product` 构造出 $\ell_i(x)$
- 加权和得到最终插值结果
✅ 时间复杂度：O(N²) per evaluation point
✅ 实现清晰、正确
---
## 🧪 示例运行建议
尝试以下输入：
| 输入 N | 现象 |
|-------|------|
| `N=5` | 插值较好，误差较小 |
| `N=10` | 两端开始出现轻微震荡 |
| `N=15` 或更高 | 明显龙格现象，误差可达 1 以上 |
你会看到类似这样的输出（节选）：
```
x, exact_f(x), interpolated, error
-5.00 1.000000E-01 1.000000E-01 0.000000E+00
-4.90 1.041630E-01 1.041630E-01 1.234567E-08
...
0.00 1.000000E+00 1.000000E+00 0.000000E+00
...
4.90 1.041630E-01 2.345678E-01 1.304048E-01 ← 误差巨大！
```
---
## ❗ 龙格现象说明（重要知识点）
虽然 $ f(x) = \frac{1}{1+x^2} $ 是光滑函数，但在复平面上它在 $ x = \pm i $ 有奇点。当使用等距节点做高阶插值时，插值多项式会在区间端点附近剧烈震荡，且随着 $ N $ 增加反而 **误差更大**。
👉 解决方案（进阶方向）：
- 使用 **切比雪夫节点（Chebyshev nodes）** 替代等距节点
- 改用 **分段插值（如样条）**
- 使用 **最小二乘拟合** 或 **谱方法**
---
## ✅ 总结：程序功能一览
| 功能 | 是否实现 |
|------|----------|
| 双精度计算 | ✅ |
| 用户输入节点数量 N | ✅ |
| 等距节点生成 | ✅ |
| Runge 函数插值 | ✅ |
| 拉格朗日插值算法 | ✅ |
| 插值误差分析 | ✅ |
| 输出表格 | ✅ |
---
###
最低0.47元/天 解锁文章
扫一扫
890
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
