转载

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

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

波动方程：

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

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

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

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

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

效果：

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

 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

欢迎大家一起探讨计算机图形学算法问题，邮箱：shushen@126.com

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

原文 http://www.cnblogs.com/shushen/p/5123964.html

正文到此结束