wave2.m

matlab算法集 matlab算法集

function [x,y,U] = wave2 (T,a,b,m,n,p,beta,f,g)
%----------------------------------------------------------------------
% Usage:       [x,y,U] = wave2 (T,a,b,m,n,p,beta,f,g)
%          
% Description: Use the explicit central difference method to solve the
%              following partial differential equation called the two-
%              dimensional wave equation:
%
%                 (d/dt)(d/dt)u(t,x,y) = beta*[(d/dx)(d/dx)u(t,x,y) +
%                                        (d/dy)(d/dy)u(t,x,y)]
%
%              The wave equation is to be solved over the region:
%
%                 0 < t <= T,
%                 0 < x < a,
%                 0 < y < b
%         
%              subject to the initial and boundary conditions:
%
%                 u(0,x,y)       = f(x,y),
%                 (d/dt)u(0,x,y) = g(x,y),
%                 u(t,x,y)       = 0
%
% Inputs:      T    = final solution time (T >= 0)
%              a    = upper boundary of x (a > 0)
%              b    = upper boundary of y (b > 0)
%              m    = number of steps in t (m >= 1)
%              n    = number of steps in x (n >= 1)
%              p    = number of steps in y (p >= 1)
%              beta = wave equation coefficient (beta > 0)
%              f    = string containing name of user-supplied 
%                     initial value function
%              g    = string containing name of user-supplied 
%                     initial derivative function
%
%                     The functions f and g are of the form:
%
%                        function z = f(x,y)
%                        function z = g(x,y)
%
%                     When f is called it must return the value u(0,x,y).
%                     When g is called it must return (d/dt)u(t,x,y)
%                     evaluated at t = 0$. The functions f and g are
%                     optional in the sense that either argument can be
%                     replaced with the empty string, ''. When this is 
%                     done, the corresponding function is assumed to be 
%                     zero which means that a user-supplied function is 
%                     NOT required in this case.
%
% Outputs:     x = n by 1 vector of x grid points
%              y = p by 1 vector of y grid points
%              U = n by p matrix containing final solution estimate 
%                  with U(i,j) = u(T,x(i),y(j)).
%
% Notes:       The central difference method is stable when T satisfies 
%              the following criterion:
%
%                 (T/m)^2 <= ((a/(n+1)^2 + (b/(p+1)^2)/(4*beta)
%
%              It is the responsibility of the user to pick T to ensure
%              that the method is stable. 
%----------------------------------------------------------------------
   
% Check arguments
   
   T = args (T,eps,T,1,'wave2');
   a = args (a,eps,a,2,'wave2');
   b = args (b,eps,b,3,'wave2');
   m = args (m,1,m,4,'wave2');
   n = args (n,1,n,5,'wave2');
   p = args (p,1,p,6,'wave2');
   beta = args (beta,eps,beta,7,'wave2');
   funf = getfun (f);
   fung = getfun (g);

% Check for stability 

   dt = T/m;
   dx = a/(n+1);
   dy = b/(p+1);
   U = zeros (n,p);
   Dt = sqrt((dx*dx + dy*dy)/(4*beta));
   if dt > Dt 
      fprintf ('\n\n------------------------------------------------------------');
      fprintf ('\nThe solution time T in the call to function wave2 is too');
      fprintf ('\nlarge for a stable solution.  It must satisfy:');
      fprintf ('\n\n  pow(T/m,2) <= (pow(a/(n+1),2) + pow(b/(p+1),2))/(4*beta)');
      fprintf ('\n\nThe maximum stable value is T = %.3f. You must decrease T',...
                m*Dt);
      fprintf ('\nor modify m, a, b, n, or p as needed to ensure stability.');
      fprintf ('\n------------------------------------------------------------\n');
      wait
      return
   end

% Compute x, y, U0 

   x = [dx : dx : n*dx]';
   y = [dy : dy : p*dy]';
   U0 = zeros (n,p);
   U1 = zeros (n,p);
   if funf
      for i = 1 : n
         for j = 1 : p
            U0(i,j) = feval(f,x(i),y(j));
         end
      end
   end

% Compute: U1(1) 

   hwbar = waitbar (0,'Solving Wave Equation: wave2'); 
   gx = beta*(dt/dx)^2;   
   gy = beta*(dt/dy)^2;
   gz = 2*(1 - gx - gy);
   if funf
      U1(1,1) = (gx*feval(f,x(2),y(1)) + gz*feval(f,x(1),y(1)) ...
                 + gy*feval(f,x(1),y(2)))/2;
      for j = 2 : p-1
         U1(1,j) = (gx*feval(f,x(2),y(j)) + gz*feval(f,x(1),y(j)) ...
                    + gy*(feval(f,x(1),y(j-1)) + ...
                    feval(f,x(1),y(j+1))))/2;
      end  
      U1(1,p) = (gx*feval(f,x(2),y(p)) + gz*feval(f,x(1),y(p)) ...
                 + gy*feval(f,x(1),y(p-1)))/2; 
   else
      U1(1,1:p) = 0;
   end
   if fung
      for j = 1 : p
         U1(1,j) = U(1,j) + dt*feval(g,x(1),y(j));
      end
   end 	

% Compute U1(i), 1 < i < n 

   
   for i = 2: n-1
      if funf
         U1(i,1) = (gx*(feval(f,x(i-1),y(1)) + feval(f,x(i+1),y(1))) ...
                    + gz*feval(f,x(i),y(1)) + gy*feval(f,x(i),y(2)))/2; 
         for j = 2 : p-1
            U1(i,j) = (gx*(feval(f,x(i-1),y(j)) + feval(f,x(i+1),y(j))) ...
                       + gz*feval(f,x(i),y(j)) + gy*(feval(f,x(i),y(j-1)) ...
                       + feval(f,x(i),y(j+1))))/2;
         end
         U1(i,p) = (gx*(feval(f,x(i-1),y(p)) + feval(f,x(i+1),y(p))) ...
                    + gz*feval(f,x(i),y(p)) + gy*feval(f,x(i),y(p-1)))/2;
      else
         U1(i,1:p) = 0;
      end
      if fung
         for j = 1 : p
            U1(i,j) = U1(i,j) + dt*feval(g,x(i),y(j));
         end
      end
   end
   
% Compute: U1(n) 

   if funf
      U1(n,1) = (gx*feval(f,x(n-1),y(1)) + gz*feval(f,x(n),y(1)) ...
                 + gy*feval(f,x(n),y(2)))/2;
      for j = 2 : p-1
         U1(n,j) = (gx*feval(f,x(n-1),y(j)) + gz*feval(f,x(n),y(j)) ...
                    + gy*(feval(f,x(n),y(j-1)) + ...
                    feval(f,x(n),y(j+1))))/2;
      end
      U1(n,p) = (gx*feval(f,x(n-1),y(p)) + gz*feval(f,x(n),y(p)) ...
                 + gy*feval(f,x(n),y(p-1)))/2; 
   else
      U1(n,1:p) = 0;
   end
   if fung
      for j = 1 : p
         U1(n,j) = U1(n,j) + dt*feval(g,x(n),y(j));
      end
   end 

   for k = 1 : m-1

      waitbar (k/m)
   
% Compute: U(1) 

      U(1,1) = gx*U1(2,1) + gz*U1(1,1) + gy*U1(1,2) - U0(1,1);
      for j = 2 : p-1
         U(1,j) = gx*U1(2,j) + gz*U1(1,j) + gy*(U1(1,j-1) + ...
                  U1(1,j+1)) - U0(1,j);
      end
      U(1,p) = gx*U1(2,p) + gz*U1(1,p) + gy*U1(1,p-1) - U0(1,p);

% Compute U(i), 1 < i < n 

      for i = 2 : n-1
         U(i,1) = gx*(U1(i-1,1) + U1(i+1,1)) + gz*U1(i,1) ...
                   + gy*U1(i,2) - U0(i,1);
         for j = 2 : p-1
            U(i,j) = gx*(U1(i-1,j) + U1(i+1,j))+ gz*U1(i,j) ...
                      + gy*(U1(i,j-1) + U1(i,j+1)) - U0(i,j);
         end
         U(i,p) = gx*(U1(i-1,p) + U1(i+1,p)) + gz*U1(i,p) ...
                   + gy*U1(i,p-1) - U0(i,p);
      end

% Compute: U(n) 

      U(n,1) = gx*U1(n-1,1) + gz*U1(n,1) + gy*U1(n,2) - U0(n,1);
      for j = 2 : p-1
         U(n,j) = gx*U1(n-1,j) + gz*U1(n,j) + gy*(U1(n,j-1) + ...
                  U1(n,j+1)) - U0(n,j);
      end
      U(n,p) = gx*U1(n-1,p) + gz*U1(n,p) + gy*U1(n,p-1) - U0(n,p);

      U0 = U1;
      U1 = U;
   end
   close (hwbar);
%----------------------------------------------------------------------