基于网格的波动方程模拟（Wave equation on mesh）附源码

　　波动方程是偏微分方程 (PDE) 里的经典方程，它在物理学中有大量应用并经常用来解释空间中的能量传播。波动方程是一个依赖时间的方程，它解释了系统状态是如何随着时间的推移而发生变化。在下面模拟波动方程时会使用会到拉普拉斯（Laplacian）算子，这是一个线性算子，具体形式在“网格形变算法”中有所解释。

　　波动方程：

其中：b为衰减系数，1/sqrt(a)为波传播速度，h为沿网格顶点法向移动的距离。

　　将波动方程离散化并整理后得到：

其中：dt为时间间隔，I为单位矩阵，L为离散Laplacian算子，h_i为待求的波高，h_i-1和h_i-2为前两个时刻的波高。

　　因此波动方程可以根据系统前两时刻的状态求解系统的当前状态。

效果：

function [Frame] = wave_equation(V, F)
    % vertex normal
    N = per_vertex_normals(V, F);
    
    % laplacian operator
    L = cotmatrix(V, F);
    
    % identical matrix
    nV = size(V, 1);
    I = speye(nV);

    % set parameters
    dt = 0.0075;
    a = 200;
    b = 0.5;
    A = (1 + b*dt)*I - a*dt^2*L;

    % figure
    figure('Position', [400, 400, 400, 320]);
    fh = trisurf(F, V(:,1), V(:,2), V(:,3), ...
        'FaceColor', 'interp', 'FaceLighting', 'phong', ...
        'EdgeColor', 'none', 'CData', zeros(nV,1));
    view(2)
    axis equal
    axis off
    camlight

    set(gca, 'position', [0 0 1 1]);
    set(fh, 'ButtonDownFcn', @ondown);

    V0 = V;
    H = zeros(nV, 1);
    H_1 = zeros(nV, 1);
    draw_t = 0;
    i = 1;
    tic;
    while true
        H_2 = H_1;
        H_1 = H;

        B = (2 + b*dt)*H_1 - H_2;
        H = AB;

        % updata figure
        draw_t = draw_t + toc;
        if draw_t > 0.033
            V = V0 + bsxfun(@times, H, N);
            set(fh, 'Vertices', V, 'CData', H);
            drawnow;

            Frame(i) = getframe(gcf);
            i = i + 1;
            if i > 150
                break;
            end
            
            tic;
            draw_t = 0;
        end
    end

    function ondown(src, ev)
        pt = get(gca, 'CurrentPoint');

        dH = 1/10*sparse(knnsearch(V, pt(1,:)), 1, 1, nV, 1);
        H = H + dH;
    end
end

本文为原创，转载请注明出处：http://www.cnblogs.com/shushen。