wave1.m

用中心差分实现一维小波变换

function [t,x,U] = wave1 (T,a,m,n,beta,f,g)
%----------------------------------------------------------------------
% Usage:       [t,x,U] = wave1 (T,a,m,n,beta,f,g);
%
% Description: Use the explicit central difference method to solve the
%              following partial differential equation called the
%              one-dimensional wave equation:
%
%                 (d/dt)(d/dt)u(t,x) = beta*(d/dx)(d/dx)u(t,x)
%
%              The wave equation is to be solved over the region:
%
%                 0 <= t <= T,
%                 0 < x < a
%
%              The boundary conditions are:
%
%                 u(t,0) = 0,
%                 u(t,a) = 0,
%
%              The initial conditions are:
%
%                 u(0,x)       = f(x)
%                 (d/dt)u(0,x) = g(x)
%
% Inputs:      T    = final solution time (T > 0)
%              a    = upper boundary of x (a > 0)
%              m    = number steps of t (m >= 1)
%              n    = number of steps of x (n >= 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 y = f(x)
%                        function y = g(x)
%
%                     When f is called it must return the initial value
%                     u(0,x). When g is called it must return (d/dt)u(t,x)
%                     evaluated at t = 0.  The functions f and g are
%                     optional in the sense that either argument can be
%                     replaced by the empty string, '', When this is
%                     done,the corresponding function is assumed to be
%                     zero which means a user-supplied function is NOT
%                     required in this case.
%
% Outputs:     t = (m+1) by 1 vector containing t grid
%              x = n by 1 vector containing x grid       
%              U = (m+1) by n matrix containing the estimated solution
%                   with U(k,j) = u(t(k),x(j).
%
% Notes:       The central difference method is stable for gamma <= 1
%              where the gain gamma is defined:
%
%                 gamma = beta*(T*(n+1)/(a*m))^2
%
%              It is the responsibility of the user to pick T, a, n,
%              and m to ensure that gamma <= 1.  
%----------------------------------------------------------------------
   
% Check arguments

   T = args (T,eps,T,1,'wave1');
   a = args (a,eps,a,2,'wave1');
   m = args (m,1,m,3,'wave1');
   n = args (n,1,n,4,'wave1');
   beta = args (beta,eps,beta,5,'wave1');
   funf = getfun (f);
   fung = getfun (g);

% Initialize 

   dt = T/m;
   dx = a/(n+1);
   t = [0 : dt : T]';
   x = [dx : dx : n*dx]'; 
   gamma =  beta*(dt/dx)^2;
   tmax =  m*dx/sqrt(beta);
   if T > tmax 
      fprintf ('\n\n--------------------------------------------------------');
      fprintf ('\nThe solution time T in the call to function wave1 is too');
      fprintf ('\nlarge for a stable solution.  It must satisfy:');
      fprintf ('\n\n    (T/(m+1)^2 <= (a/(n+1))^2/beta');
      fprintf ('\n\nThe maximum stable value is T = %.3f. You must decrease',...
                tmax);
      fprintf ('\nT or modify m, a or n as needed to ensure stability.');
      fprintf ('\n--------------------------------------------------------\n');
      wait 
      return; 
   end
   
% Initialize U(1,:) and U(2,:) using initial condition functions 

   g0 = 2*(1 - gamma);
   if funf
      for i = 1 : n
         U(1,i) = feval(f,x(i));
      end
      U(2,1)  = (g0*feval(f,x(1)) + gamma*feval(f,x(2)))/2;
      for i = 2 : n-1
         U(2,i) = (gamma*feval(f,x(i-1)) + g0*feval(f,x(i)) + ...
                   gamma*feval(f,x(i+1)))/2;
      end
      U(2,n)  = (gamma*feval(f,x(n-1)) + g0*feval(f,x(n)))/2;
   else 
      U(1:2,1:n) = 0;
   end 
   if fung
      for i = 1 : n
         U(2,i) = U(2,i) + dt*feval(g,x(i)); 
      end
   end

% Compute remaining rows of U 
             
   for i = 2 : m
      U(i+1,1)  = g0*U(i,1) + gamma*U(i,2) - U(i-1,1);
      for j = 2 : n-1
         U(i+1,j) = gamma*U(i,j-1) + g0*U(i,j) + ... 
                    gamma*U(i,j+1) - U(i-1,j);
      end 
      U(i+1,n)  = gamma*U(i,n-1) + g0*U(i,n) - U(i-1,n); 
   end
%----------------------------------------------------------------------