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
结果如下:
这个看着就挺像波动的。
和三维热传导方程类似,三维波动方程也难以画出来,这里就不再实现了。