稀疏网格高斯-埃尔米特数值积分方法

稀疏网格方法是一种高效的高维数值积分技术，通过组合一维积分规则来避免"维度灾难"。MATLAB实现用于生成稀疏网格高斯-埃尔米特积分节点和权重。

function [nodes, weights] = sparseGridGH(d, level, varargin)
% SPARSEGRIDGH 生成稀疏网格高斯-埃尔米特积分节点和权重
%   [nodes, weights] = sparseGridGH(d, level) 
%   生成d维、指定水平的稀疏网格高斯-埃尔米特积分规则
%
%   输入:
%       d      - 维度 (正整数)
%       level  - 稀疏网格水平 (正整数, 通常≥1)
%
%   输出:
%       nodes   - 积分节点 (每列是一个d维点)
%       weights - 积分权重 (行向量)
%
%   可选参数:
%       'maxPrecision' - 最大精度水平 (默认: 100)
%       'tol'         - 权重容差 (默认: 1e-10)

    % 处理可选参数
    p = inputParser;
    addParameter(p, 'maxPrecision', 100, @isnumeric);
    addParameter(p, 'tol', 1e-10, @isnumeric);
    parse(p, varargin{:});
    maxPrecision = p.Results.maxPrecision;
    tol = p.Results.tol;
    
    % 检查输入有效性
    if d < 1 || level < 1
        error('维度和水平必须为正整数');
    end
    
    % 获取一维规则序列
    oneDRules = get1DGHRules(maxPrecision);
    
    % 初始化节点和权重
    nodes = [];
    weights = [];
    
    % 计算稀疏网格
    for k = level:level+d-1
        % 生成所有可能的索引组合
        indexSets = generateIndexSets(d, k);
        
        for i = 1:size(indexSets, 1)
            indexSet = indexSets(i, :);
            
            % 检查索引集是否有效
            if any(indexSet > maxPrecision)
                warning('索引超出最大精度水平，结果可能不精确');
                continue;
            end
            
            % 计算当前索引集的组合系数
            coeff = (-1)^(k - sum(indexSet)) * nchoosek(d-1, k - sum(indexSet));
            
            % 获取当前索引集对应的规则
            rule = getProductRule(oneDRules, indexSet);
            
            % 应用组合系数并添加到网格
            if ~isempty(rule.nodes) && ~isempty(rule.weights)
                % 添加节点
                nodes = [nodes, rule.nodes];
                
                % 添加权重（应用组合系数）
                weights = [weights, coeff * rule.weights];
            end
        end
    end
    
    % 合并重复节点（考虑数值容差）
    [nodes, weights] = mergeNodes(nodes, weights, tol);
    
    % 归一化权重（确保权重和为1）
    weights = weights / sum(weights);
end

function oneDRules = get1DGHRules(maxLevel)
% GET1DGHRULES 生成一维高斯-埃尔米特积分规则序列
%   生成从1到maxLevel的一维高斯-埃尔米特积分规则
    oneDRules = cell(1, maxLevel);
    
    for l = 1:maxLevel
        n = get1DGHNumPoints(l);
        [x, w] = GaussHermite(n);
        oneDRules{l} = struct('nodes', x, 'weights', w);
    end
end

function n = get1DGHNumPoints(level)
% GET1DGHNUMPOINTS 根据水平获取一维高斯-埃尔米特积分点数
%   使用指数增长策略：点数 = 2^level - 1
    n = 2^level - 1;
    if n < 1
        n = 1;
    end
end

function [x, w] = GaussHermite(n)
% GAUSSHERMITE 生成一维高斯-埃尔米特积分点和权重
%   使用特征值方法计算n点高斯-埃尔米特积分规则
    
    % 创建对称三对角矩阵
    beta = sqrt((1:n-1)/2);
    J = diag(beta, 1) + diag(beta, -1);
    
    % 计算特征值和特征向量
    [V, D] = eig(J);
    
    % 提取节点和权重
    x = diag(D)';
    w = V(1, :).^2 * sqrt(pi);
    
    % 排序节点
    [x, idx] = sort(x);
    w = w(idx);
    
    % 调整权重符号以确保为正
    w = abs(w);
end

function indexSets = generateIndexSets(d, k)
% GENERATEINDEXSETS 生成满足 |i| = k 的索引集
%   生成所有满足 i1 + i2 + ... + id = k 的正整数索引组合
    
    % 使用递归生成所有可能的组合
    indexSets = [];
    current = ones(1, d);
    
    function generate(idx, sumLeft)
        if idx > d
            if sumLeft == 0
                indexSets = [indexSets; current];
            end
            return;
        end
        
        for val = 1:sumLeft
            current(idx) = val;
            generate(idx + 1, sumLeft - val);
        end
    end

    generate(1, k);
end

function productRule = getProductRule(oneDRules, indexSet)
% GETPRODUCTRULE 获取多维张量积规则
%   根据索引集创建多维张量积积分规则
    
    d = length(indexSet);
    rules = cell(1, d);
    
    % 获取每个维度的一维规则
    for i = 1:d
        level = indexSet(i);
        if level <= length(oneDRules)
            rules{i} = oneDRules{level};
        else
            rules{i} = struct('nodes', [], 'weights', []);
        end
    end
    
    % 初始化张量积规则
    productRule = struct('nodes', [], 'weights', []);
    
    % 递归计算张量积
    function tensorProduct(dim, currentNodes, currentWeight)
        if dim > d
            productRule.nodes = [productRule.nodes, currentNodes'];
            productRule.weights = [productRule.weights, currentWeight];
            return;
        end
        
        rule = rules{dim};
        n = length(rule.nodes);
        
        for i = 1:n
            newNodes = [currentNodes, rule.nodes(i)];
            newWeight = currentWeight * rule.weights(i);
            tensorProduct(dim + 1, newNodes, newWeight);
        end
    end

    tensorProduct(1, [], 1);
end

function [mergedNodes, mergedWeights] = mergeNodes(nodes, weights, tol)
% MERGENODES 合并重复节点（考虑数值容差）
%   合并空间中接近的节点，并累加它们的权重
    
    if isempty(nodes)
        mergedNodes = [];
        mergedWeights = [];
        return;
    end
    
    % 初始化
    n = size(nodes, 2);
    d = size(nodes, 1);
    mergedNodes = zeros(d, 0);
    mergedWeights = [];
    merged = false(1, n);
    
    % 遍历所有节点
    for i = 1:n
        if merged(i)
            continue;
        end
        
        % 查找接近的节点
        closeIndices = i;
        for j = i+1:n
            if ~merged(j) && all(abs(nodes(:, i) - nodes(:, j)) < tol)
                closeIndices = [closeIndices, j];
                merged(j) = true;
            end
        end
        
        % 合并节点和权重
        mergedNode = mean(nodes(:, closeIndices), 2);
        mergedWeight = sum(weights(closeIndices));
        
        % 添加到结果
        mergedNodes = [mergedNodes, mergedNode];
        mergedWeights = [mergedWeights, mergedWeight];
    end
end

function visualizeSparseGrid(nodes, weights, dims)
% VISUALIZESPARSEGRID 可视化稀疏网格
%   绘制二维或三维稀疏网格
    
    if nargin < 3
        dims = [1, 2];
        if size(nodes, 1) > 2
            dims = [1, 2, 3];
        end
    end
    
    % 提取要可视化的维度
    visNodes = nodes(dims, :);
    d = length(dims);
    
    % 创建图形
    figure('Color', 'white', 'Position', [100, 100, 800, 600]);
    
    if d == 2
        % 二维可视化
        scatter(visNodes(1, :), visNodes(2, :), 80 * weights / max(weights), 'filled');
        colormap(jet);
        colorbar;
        title('稀疏网格高斯-埃尔米特积分节点 (尺寸表示权重)');
        xlabel(['维度 ', num2str(dims(1))]);
        ylabel(['维度 ', num2str(dims(2))]);
        grid on;
        
    elseif d == 3
        % 三维可视化
        scatter3(visNodes(1, :), visNodes(2, :), visNodes(3, :), ...
                80 * weights / max(weights), weights, 'filled');
        colormap(jet);
        colorbar;
        title('稀疏网格高斯-埃尔米特积分节点 (尺寸和颜色表示权重)');
        xlabel(['维度 ', num2str(dims(1))]);
        ylabel(['维度 ', num2str(dims(2))]);
        zlabel(['维度 ', num2str(dims(3))]);
        grid on;
        rotate3d on;
        
    else
        warning('只能可视化2D或3D投影');
    end
end

function testSparseGridGH()
% TESTSPARSEGRIDGH 测试稀疏网格高斯-埃尔米特积分实现
    
    % 测试1: 一维情况
    fprintf('测试1: 一维高斯-埃尔米特积分\n');
    d = 1;
    level = 3;
    [nodes, weights] = sparseGridGH(d, level);
    
    % 计算简单积分 (e^{-x^2} 在 R 上的积分应为 sqrt(pi))
    f = ones(size(nodes)); % 被积函数 f(x) = 1
    integral_approx = dot(weights, f);
    integral_exact = sqrt(pi);
    fprintf('  近似积分: %.8f\n', integral_approx);
    fprintf('  精确积分: %.8f\n', integral_exact);
    fprintf('  相对误差: %.4e\n\n', abs(integral_approx - integral_exact)/integral_exact);
    
    % 测试2: 二维情况
    fprintf('测试2: 二维高斯-埃尔米特积分\n');
    d = 2;
    level = 3;
    [nodes, weights] = sparseGridGH(d, level);
    
    % 计算积分 (e^{-(x^2+y^2)} 在 R^2 上的积分应为 pi)
    f = ones(1, size(nodes, 2)); % 被积函数 f(x,y) = 1
    integral_approx = dot(weights, f);
    integral_exact = pi;
    fprintf('  近似积分: %.8f\n', integral_approx);
    fprintf('  精确积分: %.8f\n', integral_exact);
    fprintf('  相对误差: %.4e\n\n', abs(integral_approx - integral_exact)/integral_exact);
    
    % 测试3: 高维情况
    fprintf('测试3: 四维高斯-埃尔米特积分\n');
    d = 4;
    level = 3;
    [nodes, weights] = sparseGridGH(d, level);
    
    % 计算积分 (e^{-||x||^2} 在 R^d 上的积分应为 pi^{d/2})
    f = ones(1, size(nodes, 2)); % 被积函数 f(x) = 1
    integral_approx = dot(weights, f);
    integral_exact = pi^(d/2);
    fprintf('  近似积分: %.8f\n', integral_approx);
    fprintf('  精确积分: %.8f\n', integral_exact);
    fprintf('  相对误差: %.4e\n\n', abs(integral_approx - integral_exact)/integral_exact);
    
    % 可视化二维网格
    d = 2;
    level = 4;
    [nodes, weights] = sparseGridGH(d, level);
    visualizeSparseGrid(nodes, weights);
    
    % 可视化三维网格（二维投影）
    d = 3;
    level = 3;
    [nodes, weights] = sparseGridGH(d, level);
    visualizeSparseGrid(nodes, weights, [1, 2, 3]);
end

% 运行测试
testSparseGridGH();

方法说明

1. 稀疏网格方法原理

稀疏网格方法基于Smolyak算法，通过组合一维积分规则来构造高维积分规则：

A(q,d) = Σ_{q-d+1≤|i|≤q} (-1)^{q-|i|} * C(d-1, q-|i|) * (U_{i1}⊗...⊗U_{id})

其中：

• d 是维度
• q = level + d - 1 是稀疏网格水平
• U_i 是一维积分规则
• |i| = i1 + i2 + ... + id

2. 实现组成

(1) 一维高斯-埃尔米特积分规则

function [x, w] = GaussHermite(n)

使用特征值方法生成n点高斯-埃尔米特积分规则：

1. 构造对称三对角矩阵
2. 计算特征值和特征向量
3. 特征值给出节点，特征向量给出权重

(2) 索引集生成

function indexSets = generateIndexSets(d, k)

生成所有满足 i1 + i2 + ... + id = k 的正整数索引组合

(3) 张量积规则计算

function productRule = getProductRule(oneDRules, indexSet)

递归计算多维张量积积分规则

(4) 稀疏网格主函数

function [nodes, weights] = sparseGridGH(d, level, varargin)

1. 获取一维规则序列
2. 生成有效索引集
3. 计算组合系数
4. 组合张量积规则
5. 合并重复节点
6. 归一化权重

3. 测试函数

function testSparseGridGH()

验证实现正确性：

1. 一维高斯积分：\(∫e^{-x²}dx = √π\)
2. 二维高斯积分：\(∫∫e^{-(x²+y²)}dxdy = π\)
3. 四维高斯积分：\(∫e^{-||x||²}dx = π^{d/2}\)

使用

1. 生成二维稀疏网格

d = 2;       % 维度
level = 4;   % 稀疏网格水平
[nodes, weights] = sparseGridGH(d, level);

% 可视化
visualizeSparseGrid(nodes, weights);

2. 计算高维积分

d = 5;       % 维度
level = 3;   % 稀疏网格水平
[nodes, weights] = sparseGridGH(d, level);

% 定义被积函数 (示例: 高斯函数)
f = exp(-sum(nodes.^2, 1));

% 计算积分近似值
integral_approx = dot(weights, f);

3. 调整精度

% 设置更高精度
[nodes, weights] = sparseGridGH(d, level, 'maxPrecision', 200, 'tol', 1e-12);

参考代码 用于数值积分的稀疏网格高斯埃尔米特方法 www.youwenfan.com/contentcnj/84808.html

可视化说明

程序包含可视化函数，可展示二维或三维稀疏网格：

• 点的大小表示权重大小
• 颜色表示权重值（使用jet色谱）
• 支持选择特定维度进行可视化

此实现适用于高维数值积分问题，特别是在金融工程、统计物理和机器学习中的期望计算。