home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / DIFEQDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.6 KB  |  210 lines
  1. % DIFEQDEM    Demonstrate solution of ODE by collocation.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3. %                                       latest change:     June 3, '90
  4. clg;clc;home;echo on
  5. % We consider the nonlinear singularly perturbed problem
  6. %
  7. %  epsilon g''(x) + (g(x))^2 = 1  on  [0,1]
  8. %                   g'(0) = g(1) = 0 .
  9. %
  10.  
  11.    pause                % Touch any key to continue
  12.  
  13. % We seek an approximate solution by collocation from C^1 piecewise
  14. % cubics with a suitable breakpoint sequence; e.g.,
  15.  
  16. breaks = [0:4]/4;
  17.  
  18. % hence with the knot sequence
  19.  
  20. knots = sort([0 0 breaks breaks 1 1])
  21.  
  22. % which has the same effect as the statement 
  23. %    knots = augknt(breaks,4,2)
  24.    pause                % Touch any key to continue
  25. echo off       % for a quick check
  26. zcheck(knots-augknt(breaks,4,2));echo on
  27.  
  28.  
  29. %     This gives a quadruple knot at both 0 and 1, which is consistent with
  30. % the fact that we have cubics, i.e., have order
  31.  
  32. k = 4;
  33.  
  34. % This implies that we have
  35.  
  36. n = length(knots)-k
  37.  
  38. % degrees of freedom. This fits nicely with the fact that we expect to 
  39. % collocate at two points per polynomial piece, for a total of 8 conditions,
  40. % bringing us to 10 conditions altogether once we add the two side conditions.
  41.  
  42.    pause                % Touch any key to continue
  43.  
  44. %    We choose two Gauss points for each interval. For the `standard'
  45. % interval [-1,1]/2 of unit length, these are the two points
  46.  
  47. gauss = .5773502692*[-1 1]/2;
  48.  
  49. %    From this, we obtain the whole collection of collocation points by
  50.  
  51. ninterv = length(breaks)-1;
  52. (breaks(2:ninterv+1)+breaks(1:ninterv))/2;
  53. ans'*[1 1] + diff(breaks)'*gauss;
  54. colpnts = sort(ans(:)');
  55.  
  56.    pause                % Touch any key to continue
  57.  
  58. %    The numerical problem we want to solve is to find a pp  y  of the given
  59. % order and with the given knots which satisfies the (nonlinear) system
  60. %     
  61. %                           y'(0)  =  0
  62. %       (y(x))^2 + epsilon y''(x)  =  1   for  x in  colpnts
  63. %                            y(1)  =  0
  64.  
  65. % If  y  is our current approximation to the solution, then the linear
  66. % problem for the better (?) solution  z  by Newton's method reads
  67. %     
  68. %                           z'(0)  =  0
  69. %     w_0(x)z(x) + epsilon z''(x)  =  b(x)   for  x in  colpnts
  70. %                            z(1)  =  0
  71. %
  72. % with  w_0(x) : =  2y(x), b(x) : =  (y(x))^2 + 1 .
  73.  
  74.    pause                % Touch any key to continue
  75.  
  76. % In fact, by choosing   w_0(1): = 1, w_1(0): = 1 , and
  77. %
  78. %     w_2(x) : =  epsilon,    w_1(x) = 0    for  x  in  colpnts
  79. %
  80. % and choosing all other values of  w_0, w_1, w_2, b  not yet specified
  81. % to be zero, we can give our system the uniform shape
  82.  
  83. %  w_0(x)z(x) + w_1(x)z'(x) + w_2(x)z''(x) = b(x)   for  x  in  points
  84.  
  85. % with
  86.  
  87. points = [0,colpnts,1];
  88.  
  89.    pause                % Touch any key to continue
  90.  
  91. %   This system converts into one for the B-spline coefficients of its
  92. % solution  z . For this, we need the zeroth, first, and second derivative
  93. % at every  x  in  points  and for every relevant B-spline. These values
  94. % are supplied by the M-file  spcol .
  95.  
  96. help spcol
  97.  
  98.    pause                % Touch any key to continue
  99.  
  100. % We use  spcol  to supply the matrix
  101. colmat = spcol(knots,k,sort([points points points]));
  102.  
  103. % From this, we get the collocation matrix by combining the row triple 
  104. % associated with  x  using the weights  w_0(x), w_1(x), w_2(x)  to get the 
  105. % row corresponding to  x  of the matrix of our linear system. 
  106.  
  107.    pause                % Touch any key to continue
  108.  
  109.  
  110. % For this, we need to choose  epsilon . The problem is already very 
  111. % difficult for  epsilon = .001 , so we choose the mild
  112.  
  113. epsilon = .1;
  114.  
  115. %     We also need a current approximation  y . Initially, we get it
  116. % by interpolating some reasonable initial guess from our pp space at all
  117. % the  points .  We use the function
  118. %   ()^2 - 1
  119. % which does satisfy the end conditions, and pick the matrix from the full
  120. % matrix  colmat . Here it is, in several (cautious) steps:
  121. intmat = colmat([2 1+[1:(n-2)]*3,1+(n-1)*3],:);
  122. coefs = intmat\[0 colpnts.*colpnts-1 0]';
  123. y = spmak(knots,coefs');
  124. % We plot the result (it should be exactly  ()^2-1 ), to be sure:
  125. fnplt(y);pause
  126. hold on
  127. % That's ok, then.
  128.  
  129.    pause                % Touch any key to continue
  130.  
  131. %    We can now complete the construction and solution of the linear 
  132. % system for the improved approximate solution  z  from our current 
  133. % guess  y . In fact, with the initial guess  y  available, we now set 
  134. % up an iteration, to be terminated when the change  z-y  is `small
  135. % enough' .  We choose the tolerance
  136.  
  137. tolerance = 6.e-9;
  138.  
  139.    pause                % Touch any key to continue
  140.  
  141. difeqite
  142.  
  143. % that looks like quadratic convergence, as expected.
  144.  
  145.    pause                % Touch any key to continue
  146.  
  147. echo off %  (for a check of newknt etc. that will take a moment)
  148. typ =[-0.97210037583603;
  149.       -0.97138283316330;
  150.       -0.96077877735438;
  151.       -0.94305590383320;
  152.       -0.89729018696592;
  153.       -0.83892108187956;
  154.       -0.70191870276411;
  155.       -0.53656329907016;
  156.       -0.17862024561499;
  157.        0];
  158. zcheck(typ'-fnval(y,points));
  159.  
  160.    difeqset
  161.  
  162. echo off       % for a quick check
  163.    typ=[-0.97209501816518;
  164.         -0.97162227476240;
  165.         -0.96470968097732;
  166.         -0.95393233409872;
  167.         -0.93022030933061;
  168.         -0.90278210186314;
  169.         -0.84585242261528;
  170.         -0.78262837036445;
  171.         -0.64997520321340;
  172.         -0.49555677657819;
  173.         -0.16116850509888;
  174.                         0];
  175.    zcheck(fnval(y,points)'-typ);
  176. echo on
  177.  
  178. %  If we now decrease  epsilon , we create more of a boundary layer near 
  179. % the right endpoint, and this calls for a nonuniform mesh. We use  newknt
  180. % to construct an appropriate (finer) mesh from the current approximation, 
  181. % then divide  epsilon  by 3 and repeat the above calculation.
  182.  
  183. while 1,
  184.    echo on
  185.    difeqset  % generates new breaks, new collocation points, and
  186.              % constructs a new approximate solution by interpolating the
  187.              % current approximation from the new pp space.
  188.  
  189.    % We plot the result (it should be close to our current solution) to be sure:
  190.    fnplt(y,'x');pause
  191.    epsilon = epsilon/3
  192.    difeqite  % runs the Newton iteration 
  193.  
  194.    echo off
  195.    fprintf(' If you want to decrease  epsilon  further,\n')
  196.    fprintf(' type any digit, then RETURN. ... (then wait) ...')
  197.      answer = input(' Otherwise, just hit RETURN :  ');
  198.    if isempty(answer), break,end
  199. end
  200. echo on
  201.    
  202. % Here is the final distribution of breaks
  203.  
  204. breaks = ppbrk(sp2pp(z));
  205. for j = 1:length(breaks)
  206.    plot(breaks(j)*[1 1],[0,-1])
  207. end
  208.  
  209. hold off
  210.