基於網格的波動方程模擬（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 = A\B;

        % 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。