wave2.m

matlab一些典型例程

function [u,x,y,t] = Wave2(A,D,T,it0,i1t0,bxyt,Mx,My,N)
%解方程 a(u_xx + u_yy) = u_tt for D(1) <= x <= D(2), D(3) <= y <= D(4), 0 <= t <= T
% 初始条件: u(x,y,0) = it0(x,y), u_t(x,y,0) = i1t0(x,y)
% 边界条件: u(x,y,t) = bxyt(x,y,t) for (x,y) on Boundary
% Mx/My  x轴和y轴的等分段数
% N = 时间轴的等分段数
dx = (D(2)-D(1))/Mx; x = D(1)+[0:Mx]*dx;
dy = (D(4)- D(3))/My; y = D(3)+[0:My]'*dy;
dt = T/N; t = [0:N]*dt;
%初始化
u = zeros(Mx+1,My + 1); ut = zeros(Mx+ 1,My + 1);
for i = 2:Mx
    for j = 2:My
        u(i,j) = it0(x(i),y(j)); ut(i,j) = i1t0(x(i),y(j));
    end
end
adt2 = A*dt*dt; rx = adt2/(dx*dx); ry = adt2/(dy*dy);
rxy1 = 1- rx - ry; rxy2 = rxy1*2;
u_1 = u;
for k = 0:N
    t = k*dt;
    for i = 1:Mx + 1 %边界条件
        u(i,[1 My + 1]) = [bxyt(x(i),y(1),t) bxyt(x(i),y(My+1),t)];
    end
    for j = 1:My + 1
        u([1 Mx + 1],j) = [bxyt(x(1),y(j),t); bxyt(x(Mx+1),y(j),t)];
    end
    if k == 0
        for i = 2:Mx
            for j = 2:My %(11.3.8)
                u(i,j) = 0.5*(rx*(u_1(i-1,j) + u_1(i+1,j))+ ry*(u_1(i,j-1)+u_1(i,j-1))) + rxy1*u(i,j) + dt*ut(i,j);
            end
        end
    else
        for i = 2:Mx
            for j = 2:My %11.3.7
                u(i,j) = rx*(u_1(i-1,j)+ u_1(i+1,j))+ ry*(u_1(i,j-1) + u_1(i,j+1)) + rxy2*u(i,j) -u_2(i,j);
            end
        end
    end
    u_2 = u_1; u_1 = u;
end
u=u';