Math 563 Scientific Computing II

Homework 18 Due Mar.26 Hangjie Ji

Question: Program the deferred correction algorithm to solve the homework problem for March 7. $y'(t) = f(y(t))$, where y and f are 2-vectors. $f_0 = 2 y_0 ( 1 - y_1 )$ and $f_1 = - y_1 ( 1 - y_0 )$. The initial condition is $y_0(0) = 1$ and $y_1(0) = 3$. Use a second-order Runge-Kutta scheme (such as the modified Euler method, improved Euler method, or the SDIRK scheme) to solve the ordinary differential equations. Choose parameters in the deferred correction scheme to achieve order 12, and describe how you chose those parameters. Also describe how you chose your timestep for the deferred correction method. Perform a mesh refinement study to verify that your method achieves the desired order. You may use the code in deferred_correction.f and GUIDeferredCorrection.C to help you.

Discussion:

I used the modified Euler method to solve the ordinary differential equations which looks like \[ y_{i+\frac{1}{2}} = y_i + \frac{t_{i+1} - t_{i}}{2}f(y_i) \] \[ y_{i+1} = y_i + (t_{i+1} - t_{i})f(y_{i+\frac{1}{2}}). \] To achieve order 12, in each time step the spectral deferred correction should consist of one second order prediction (modified Euler method, for instance) and five subsequent deferred corrections, each of which takes the form \[ \epsilon = \int_0^{t_i} L^m(\{ f(y_i, t_i) \},\tau)d\tau - y_{i+1} + y_0 \] \[ d_{i+\frac{1}{2}} = d_i + \frac{t_{i+1} - t_{i}}{2}(f(y_i + d_i,t_i) - f(y_i,t_i)) \] \[ d_{i+1} = d_i + (t_{i+1} - t_{i})(f(d_{i+\frac{1}{2}}+d_i,t_i) - f(d_{i+\frac{1}{2}},t_i)) + [\epsilon_{i+1} - \epsilon_{i}]. \] For the nonlinear initial value problem \[ y_0' = 2 y_0 ( 1 - y_1 ) \] \[ y_1' = - y_1 ( 1 - y_0 ) \] \[ y_0(0) = 1 , y_1(0) = 3, \] The corresponding Jacobian matrix is

\begin{bmatrix} 2-2y_1 & -2y_0\\ y_1 & -1+y_0. \end{bmatrix}

I used fixed time step in this problem. But generally to choose an appropriate timestep for the deferred correction method, the product of the timestep and the maximum decay mode should lie in the region of absolute stability of the method which is used as a prediction.

I find it hard to compile the mentioned codes in C and Fortran, so I simply use Matlab to implement the discussed deferred correction method. The Lagrange interpolating polynomials are computed using five points Gaussian quadratures. The Lobatto quadrature points ${s_i}$ on the interval [-1,1], as wells as Gaussian quadrature points and corresponding weights are defined directly in the code for convenience. Finally the mesh refinement study (Figure 2) shows that the convergence is order 12.

Output:

Figure 1

Figure 2

The source code is attached here: C++ code colored by C++2HTML 
  close all;
  clear all;
  clc;
 
   
  % Lobatto quadrature points on [-1,1];
  s = [-1,-sqrt((7+2*sqrt(7))/21),-sqrt((7-2*sqrt(7))/21),sqrt((7-2*sqrt(7))/21),sqrt((7+2*sqrt(7))/21),1]; 
  
  In = zeros(5,6);
  weight = [(322-13*sqrt(70))/900,(322+13*sqrt(70))/900,128/225,(322+13*sqrt(70))/900,(322-13*sqrt(70))/900];
  a = [-sqrt(5+2*sqrt(10/7))/3,-sqrt(5-2*sqrt(10/7))/3,0,sqrt(5-2*sqrt(10/7))/3,sqrt(5+2*sqrt(10/7))/3];
  
  for k = 1:5
      q = (s(k+1)-s(1))/2*a + (s(k+1)+s(1))/2;
      
      for i = 1:6
          % compute for In(k,i)
          for i_quad = 1 : 5     
              f = 1;
               for j = 1 : 6
                    if j == i
                        continue;
                    end
                    f = f *(q(i_quad) - s(j))/(s(i) - s(j));
                    
               end
               In(k,i) = In(k,i) + weight(i_quad) * f;
          end
          In(k,i) = In(k,i)*(s(k+1)-s(1))/2;
      end
  end
  
  % mesh refinement study
  y_final = zeros(7,2);
  h = 0.8; % timestep
  h_timestep = zeros(7,1);
  for n_iteration = 1:7
      h_timestep(n_iteration) = h;
      
      step_num = floor(10/h);
      
      t = h*0.5 + h*0.5*s;

      y_old = [1,3];
      y = zeros(6,2);
      y_record = zeros(step_num,2);

      for t_step = 1 : step_num
         
          y(1,:) = y_old;
          for i = 1 : 5

              tmp = y(i,:) + 0.5*(t(i+1) - t(i))*objfun(y(i,:));
              y(i+1,:) = y(i,:) + (t(i+1) - t(i))*objfun(tmp);
          end

          d = zeros(6,2);
          err = zeros(6,2);
          for j = 1 : 6
              err(1,:) = [0,0];
              for i = 1 : 5
                  err_y = [0,0];

                  for m = 1 : 6
                      err_y = err_y + In(i,m)*objfun(y(m,:));
                  end

                  err(i+1,:) = 0.5*h*err_y - y(i+1,:) + y(1,:);
                  
                  tmp2 = d(i,:) + 0.5*(t(i+1) - t(i))* (objfun(y(i,:)+d(i,:))- objfun(y(i,:)));
                  d(i+1,:) = d(i,:) + (t(i+1) - t(i))* (objfun(tmp2 + d(i,:)) - objfun(tmp2)) + err(i+1,:) - err(i,:);

              end

              for k = 2 : 6

                  y(k,:) = y(k,:) + d(k,:);
              end
          end
          
          y_old = y(end,:);
           y_record(t_step,:) = y_old;

      end
         
      y_final(n_iteration,:) = y_record(end,:);
      h = h/2; % mesh refinement
  end
  
 
  
  figure(1),plot(y_record(:,1),y_record(:,2)),xlabel('y1','fontsize',14), 
  ylabel('y2','fontsize',14),title('Phase Portrait','fontsize',14);
 
  y_exact = y_final(end,:);
  error(:,1) = abs(y_final(1:end-1,1) - y_exact(1));
  error(:,2) = abs(y_final(1:end-1,2) - y_exact(2));
  figure(2),plot(log10(h_timestep(1:end-1)),log10(error(:,1)),...
      log10(h_timestep(1:end-1)),log10(error(:,2))),xlabel('log10(timestep)','fontsize',14),
  ylabel('log10(error)','fontsize',14), title('mesh refinement study','fontsize',14);

Last modified: