matlab練習程序（差分法解二維波動方程）

上一篇實現了一維波動方程數值解，這一篇實現二維波動方程數值解。

二維波動方程如下：

寫成差分形式：

 

整理一下就能得到u(i+1,j,k)。

matlab代碼如下：

clear all;close all;clc;

t = 3;          %時間范圍，計算到3秒
x = 1;y = 1;    %空間范圍，0-1米
m = 320;        %時間t方向分320個格子
n = 32;         %空間x方向分32個格子
k = 32;         %空間y方向分32個格子
ht = t/(m-1);   %時間步長dt
hx = x/(n-1);   %空間步長dx
hy = y/(k-1);   %空間步長dy

u = zeros(m,n,k);

%設置邊界
[x,y] = meshgrid(0:hx:1,0:hy:1);
u(1,:,:) = sin(4*pi*x)+cos(4*pi*y);
u(2,:,:) = sin(4*pi*x)+cos(4*pi*y);

%按照公式進行差分
for ii=2:m-1
    for jj=2:n-1
        for kk=2:k-1
            u(ii+1,jj,kk) = ht^2*(u(ii,jj+1,kk)+u(ii,jj-1,kk)-2*u(ii,jj,kk))/hx^2 + ...
                ht^2*(u(ii,jj,kk+1)+u(ii,jj,kk-1)-2*u(ii,jj,kk))/hy^2 + 2*u(ii,jj,kk) - u(ii-1,jj,kk);
        end
    end
end

for i=1:320
    figure(1);
    mesh(x,y,reshape(u(i,:,:),[n k]));
    axis([0 1 0 1 -4 4]);
    
%     F=getframe(gcf);
%     I=frame2im(F);
%     [I,map]=rgb2ind(I,256); 
%     if i == 1
%         imwrite(I,map,'test.gif','gif','Loopcount',inf,'DelayTime',0.03);
%     else
%         imwrite(I,map,'test.gif','gif','WriteMode','append','DelayTime',0.03);    
%     end
end

結果如下：

這個看着就挺像波動的。

和三維熱傳導方程類似，三維波動方程也難以畫出來，這里就不再實現了。