引入限制器MUSCL算法示例

MUSCL算法是一种高分辨率有限体积方法，广泛应用于计算流体力学、天体物理学、弹性力学等领域。该算法的核心思想是通过引入限制器来减小数值解的震荡，从而提高数值解的精度和稳定性。

下面是一个简单的一维线性对流方程的MUSCL算法示例：

#include <iostream>
#include <cmath>

using namespace std;

const int N = 100;
const double dx = 0.1;
const double dt = 0.01;
const double u = 1.0;
const double cfl = 0.8;

double f(double x)  // 定义初始条件
{
    if (x >= 0.0 && x <= 1.0) return 1.0;
    else return 0.0;
}

double minmod(double a, double b)  // 定义 minmod 函数
{
    if (a > 0 && b > 0) return min(a, b);
    else if (a < 0 && b < 0) return max(a, b);
    else return 0.0;
}

int main()
{
    double u[N], u_star[N], u_new[N];
    double a[N], b[N], c[N];
    double x[N], t = 0.0;
    
    // 初始化
    for (int i = 0; i < N; i++) {
        x[i] = i * dx;
        u[i] = f(x[i]);
        u_star[i] = u[i];
        u_new[i] = u[i];
    }
    
    while (t < 1.0) {
        // 计算 a,b,c
        for (int i = 1; i < N - 1; i++) {
            double ul = u[i - 1];
            double ur = u[i];
            double ux = (ur - ul) / dx;
            a[i] = minmod(2.0 * ux, ul - ur);
            ul = u[i];
            ur = u[i + 1];
            ux = (ur - ul) / dx;
            b[i] = minmod(2.0 * ux, ur - ul);
            c[i] = (u[i + 1] - u[i - 1]) / (2.0 * dx);
        }
        
        // 计算 u_star
        for (int i = 1; i < N - 1; i++) {
            double u_l = u[i] - 0.5 * minmod(a[i], c[i]);
            double u_r = u[i] + 0.5 * minmod(b[i], c[i]);
            u_star[i] = u[i] - 0.5 * (u_r - u_l) * dt / dx;
        }
        
        // 计算 u_new
        for (int i = 1; i < N - 1; i++) {
            double u_l = u_star[i - 1] + 0.5 * minmod(a[i - 1], c[i - 1]);
            double u_r = u_star[i] - 0.5 * minmod(b[i], c[i]);
            u_new[i] = u_star[i] - 0.5 * (u_r - u_l) * dt / dx;
        }
        
        // 更新 u
        for (int i = 1; i < N - 1; i++) {
            u[i] = u_new[i];
        }
        
        // 更新时间
        t += dt;
        
        // 检查 CFL 条件
        double max_speed = 0.0;
        for (int i = 0; i < N; i++) {
            max_speed = max(max_speed, abs(u[i]));
        }
        if (max_speed * dt / dx > cfl) {
            cout << "CFL 条件不满足！" << endl;
            break;
        }
    }
    
    // 输出结果
    for (int i = 0; i < N; i++) {
        cout << x[i] << " " << u[i] << endl;
    }
    
    return 0;
}

在上面的代码中，我们首先定义了初始条件函数 `f(x)`，然后初始化了网格、时间步长、速度等参数。在主循环中，我们首先计算了限制器的斜率项 `a[i]`、`b[i]` 和 `c[i]`，然后利用 MUSCL 方法计算了预测值 `u_star[i]`，最后利用 Lax-Wendroff 方法计算了修正值 `u_new[i]`，并更新了数值解 `u[i]`。在每一步更新后，我们还需要检查 CFL 条件，以保证数值解的稳定性。

需要注意的是，上述代码只是一个简单的 MUSCL 算法示例，实际应用中可能需要进行更复杂的限制器设计和参数调整。