Math 563 Scientific Computing II

Homework 15 Due Mar.7 Hangjie Ji

Question:

There are a number of test problems for initial value problem solvers at John Burkardt web page Problem 6 on this web page is y'(t) = f(y(t)), where y and f are 2-vectors. In this case, 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. Program Euler's method and the second-order Adams-Bashforth method for this problem. How should you choose the timestep size for each? Plot y_1(t) versus y_0(t) for t in [0,10] for different timestep sizes, and discuss how accurate each of these two methods is.

Discussion:

The scheme for Euler method is \[y_{n+1} = y_n + hf(y_n). \] The one for 2nd order Adams-Bashforth method looks like \[y_{n+2} = y_n+1 + h(\frac{3}{2}f(t_{n+1},y_{n+1})-\frac{1}{2}f(t_n,y_n)).\]

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}.

with eigenvalues $\lambda$ of real part -2 when $(y_0,y_1) = (1,3).$ According to Euler method absolute stability regions, we need $-2 \le \lambda h \le 0.$ According to 2nd order Adams-Bashforth absolute stability regions, we need $-1 \le \lambda h \le 0.$ Thus we may choose the size of time step h according to the value of $\lambda.$ However it's hard to bound the eigenvalues of the Jacobian matrix when y varies. At the initial point, h for the Euler method should satisfy $h \le 1$, h for 2nd order Adams-Bashforth needs to satisfy $h \le 0.5.$ But smaller time step will be required in the course of time. Besides the time step for 2nd order Adams-Bashforth should be smaller than the one for Euler method in order to attain stability. Practically I chose number of steps in the initial value problem solver to be greater than 100 (h < 0.01) for both methods.

I assume the exact solution to the problem (Figure 1) to be the numerical solution of the Matlab solver ode45 with relative tolerance 1e-8. The numerical result using both methods with different time steps are presented in Figure 2 - 4. Mesh Refinement Study is also conducted for both methods, shown in Figure 5. The accuracy of Euler Method is nearly $O(h)$, while the accuracy of the second order Adams Bashforth method is $O(h^2)$.

Output:

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

C++ code colored by C++2HTML 
EULER_TEST
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.100000

      0.000000      1.000000      3.000000
      1.000000      0.038552      1.394301
      2.000000      0.044457      0.505660
      3.000000      0.154097      0.192398
      4.000000      0.755601      0.097499
      5.000000      3.882193      0.200284
      6.000000      0.249087      5.302580
      7.000000      0.000060      1.909950
      8.000000      0.000035      0.665987
      9.000000      0.000103      0.232230
     10.000000      0.000495      0.080994

  Expected final conditions:

     10.000000      3.144351      0.348822

Adam_Bashforth_2_test
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.100000

      0.000000      1.000000      3.000000
      1.000000      0.087021      1.457880
      2.000000      0.096576      0.582529
      3.000000      0.322362      0.256896
      4.000000      1.545756      0.203877
      5.000000      4.027862      1.607727
      6.000000      0.169218      2.144245
      7.000000      0.077514      0.868940
      8.000000      0.180975      0.358980
      9.000000      0.785705      0.197626
     10.000000      3.519399      0.457758

  Expected final conditions:

     10.000000      3.144351      0.348822

EULER_TEST
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.040000

      0.000000      1.000000      3.000000
      1.000000      0.061546      1.437806
      2.000000      0.068873      0.549361
      3.000000      0.237222      0.225736
      4.000000      1.180478      0.144933
      5.000000      4.887766      0.761213
      6.000000      0.126499      2.941456
      7.000000      0.017804      1.108244
      8.000000      0.031220      0.408280
      9.000000      0.129884      0.157389
     10.000000      0.723914      0.079774

  Expected final conditions:

     10.000000      3.144351      0.348822

Adam_Bashforth_2_test
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.040000

      0.000000      1.000000      3.000000
      1.000000      0.078809      1.463465
      2.000000      0.086699      0.578661
      3.000000      0.295551      0.250395
      4.000000      1.461251      0.189595
      5.000000      4.055641      1.464814
      6.000000      0.174013      2.241086
      7.000000      0.066922      0.902463
      8.000000      0.151714      0.365335
      9.000000      0.668690      0.188728
     10.000000      3.200708      0.361607

  Expected final conditions:

     10.000000      3.144351      0.348822

EULER_TEST
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.010000

      0.000000      1.000000      3.000000
      1.000000      0.073382      1.457952
      2.000000      0.080985      0.570856
      3.000000      0.277694      0.243291
      4.000000      1.382558      0.175752
      5.000000      4.354100      1.270715
      6.000000      0.161136      2.393468
      7.000000      0.050091      0.945528
      8.000000      0.107976      0.370884
      9.000000      0.475606      0.172721
     10.000000      2.496498      0.213900

  Expected final conditions:

     10.000000      3.144351      0.348822

Adam_Bashforth_2_test
  Problem number = 6
  Problem 6, Enright and Pryce #B1
  The system is autonomous.
  Number of equations is 2
  Stepsize H = 0.010000

      0.000000      1.000000      3.000000
      1.000000      0.077433      1.464390
      2.000000      0.085081      0.577996
      3.000000      0.291168      0.249322
      4.000000      1.447456      0.187358
      5.000000      4.051935      1.441014
      6.000000      0.175507      2.257557
      7.000000      0.065401      0.908417
      8.000000      0.147479      0.366631
      9.000000      0.651607      0.187635
     10.000000      3.147579      0.349496

  Expected final conditions:

     10.000000      3.144351      0.348822

I attached the Matlab code here: 

  close all;

  main.m
  
  step_num = 1000;

  [t_e,y_0_e,y_1_e] = euler_test ( 6, step_num );
  [t,y_0,y_1] = adam_bashforth_2_test (6, step_num );

  figure(1),test_ode_test04;

  figure(2),
  plot ( t_e, y_0_e, 'b-', 'Linewidth', 3 ),hold on;
  plot ( t_e, y_1_e, 'g-', 'Linewidth', 3 );
  title ('Euler Method, 1000 step' );
    grid on
    xlabel ( '<--- T --->' );
    ylabel ( '<--- Y(T) --->' );

  figure(3),
  plot ( t, y_0, 'b-', 'Linewidth', 3 ),hold on;
  plot ( t, y_1, 'g-', 'Linewidth', 3 );
  title ('2nd order Adam Bashforth, 1000 step' );
    grid on
    xlabel ( '<--- T --->' );
    ylabel ( '<--- Y(T) --->' );

    step_num = 2^5;
    error_e = zeros(7,2);
    error_a = zeros(7,2);
    step = zeros(7,1);
    %y_0_exact = y_0(end);
    %y_1_exact = y_1(end);

    y_0_exact =  3.144351304416428;
    y_1_exact = 0.348822692323950;

   for i = 1:10
      [~,y_0_e,y_1_e] = euler_test ( 6, step_num );
      [~,y_0,y_1] = adam_bashforth_2_test (6, step_num );
      error_e(i,1) = y_0_e(end) - y_0_exact;
      error_e(i,2) = y_1_e(end) - y_1_exact;
      error_a(i,1) = y_0(end) - y_0_exact;
      error_a(i,2) = y_1(end) - y_1_exact;
      step(i) = step_num;
      step_num = step_num*2;
   end

   figure(4),
   plot(log10(step),log10(abs(error_e(:,1))),'b-', 'Linewidth', 3 ),hold on;
   plot(log10(step),log10(abs(error_e(:,2))),'g-', 'Linewidth', 3 );
    title ('Euler Method, error' );
    grid on
    xlabel ( '<--- log10(nSteps) --->' );
    ylabel ( '<--- log10(error(end)) --->' );

   figure(5),
   plot(log10(step),log10(abs(error_a(:,1))),'b-', 'Linewidth', 3 ),hold on;
   plot(log10(step),log10(abs(error_a(:,2))),'g-', 'Linewidth', 3 )
   title ('2nd Order Adam Bashforth, error' );
    grid on
    xlabel ( '<--- log10(nSteps) --->' );
    ylabel ( '<--- log10(error(end)) --->' );

function y2 = p00_adam_bashforth_2_step ( test, neqn, t0, y0, t1, y1, t2 )

%*****************************************************************************80
%
%
%  Parameters:
%
%    Input, integer TEST, the problem number.
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T0, Y0(NEQN), the arguments of the derivative function.
%
%    Input, real T1, the point at which an estimate of the solution
%    is desired.
%
%    Output, real Y1(NEQN), the estimated solution at T1.
%
  dt = t1 - t0;

  yp0 = p00_fun ( test, neqn, t0, y0 );
  yp1 = p00_fun ( test, neqn, t1, y1 );

  y2(1:neqn,1) = y1(1:neqn,1) + dt * (1.5*yp1(1:neqn,1)-0.5*yp0(1:neqn,1));

  return
end


Published with MATLAB® 7.14

function [t,y_0,y_1] = adam_bashforth_2_test ( test, step_num )


%
% Adam_Bashforth_2_test uses 2nd order Adam_Bashforth Method on a test problem.
%

  if ( nargin < 1 )
    test = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Since the user did not choose a test,\n' );
    fprintf ( 1, '  the default test %d was chosen.\n', test );
  end

  if ( nargin < 2 )
    step_num = 25;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Since the user did not choose a number of steps,\n' );
    fprintf ( 1, '  the default value %d was chosen.\n', step_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Adam_Bashforth_2_test\n' );
  fprintf ( 1, '  Problem number = %d\n', test );

  title = p00_title ( test );

  fprintf ( 1, '  %s\n', title );
%
%  Initialize any problem parameters.
%
  value = p00_param_default ( test );
  if ( 0 < value )
    fprintf ( 1, '  This problem has %d parameters.\n', value );
  end
%
%  Autonomous?
%
  if ( p00_autonomous ( test ) )
    fprintf ( 1, '  The system is autonomous.\n' );
  else
    fprintf ( 1, '  The system is not autonomous.\n' );
  end
%
%  Get the number of equations.
%
  neqn = p00_neqn ( test );

  fprintf ( 1, '  Number of equations is %d\n', neqn );

  if ( 8 < neqn )
    fprintf ( 1, '  The system is large.\n' );
    fprintf ( 1, '  Print only MIN, AVERAGE, MAX and L2NORM.\n' );
  end
%
%  Get the starting point.
%
  [ t_start, y_start ] = p00_start ( test, neqn );
%
%  Get the stopping point.
%
  [ t_stop, y_stop ] = p00_stop ( test, neqn );
%
%  Print the stepsize;
%
  fprintf ( 1, '  Stepsize H = %f\n', ( t_stop - t_start ) / step_num );

  fprintf ( 1, '\n' );
  if ( neqn <= 4 )
    fprintf ( 1, '  %12f', t_start );
    for i = 1 : neqn
      fprintf ( 1, '  %12f', y_start(i) );
    end
    fprintf ( 1, '\n' );
  elseif ( neqn <= 8 )
    fprintf ( 1, '  %12f', t_start );
    for i = 1 : 4
      fprintf ( 1, '  %12f', y_start(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '              ' );
    for i = 5 : neqn
      fprintf ( 1, '  %12f', y_start(i) );
    end
    fprintf ( 1, '\n' );
  else
    y_min = min ( y_start(1:neqn) );
    y_ave = sum ( y_start(1:neqn) ) / neqn;
    y_max = max ( y_start(1:neqn) );
    y_norm = sqrt ( sum ( y_start(1:neqn).^2 ) );
    fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', ...
      t_start, y_min, y_ave, y_max, y_norm );
  end

  y0(1:neqn,1) = y_start(1:neqn,1);
  y1(1:neqn,1) = y_start(1:neqn,1);

  % store y(t) for plot
  t = zeros(step_num,1);
  y_0 = zeros(step_num,1);
  y_1 = zeros(step_num,1);
  y_0(1) = y0(1,1);
  y_1(1) = y0(2,1);

  for step = 1 : step_num

    t0 = ( ( step_num - step + 2 ) * t_start ...
         + (            step - 2 ) * t_stop ) ...
         / ( step_num );

    t1 = ( ( step_num - step + 1 ) * t_start ...
         + (            step - 1 ) * t_stop ) ...
         / ( step_num );

    t2 = ( ( step_num - step ) * t_start ...
         + (            step ) * t_stop ) ...
         / ( step_num );

    y2 = p00_adam_bashforth_2_step ( test, neqn, t0, y0, t1,y1,t2);

    if ( mod ( 10 * step, step_num ) == 0 || step == step_num )
      if ( neqn <= 4 )
        fprintf ( 1, '  %12f', t2 );
        for i = 1 : neqn
          fprintf ( 1, '  %12f', y2(i) );
        end
        fprintf ( 1, '\n' );
      elseif ( neqn <= 8 )
        fprintf ( 1, '  %12f', t2 );
        for i = 1 : 4
          fprintf ( 1, '  %12f', y2(i) );
        end
        fprintf ( 1, '\n' );
        fprintf ( 1, '              ' );
        for i = 5 : neqn
          fprintf ( 1, '  %12f', y2(i) );
        end
        fprintf ( 1, '\n' );
      else
        y_min = min ( y2(1:neqn) );
        y_ave = sum ( y2(1:neqn) ) / neqn;
        y_max = max ( y2(1:neqn) );
        y_norm = sqrt ( sum ( y2(1:neqn).^2 ) );
        fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', ...
          t1, y_min, y_ave, y_max, y_norm );
      end
    end

    y0(1:neqn,1) = y1(1:neqn,1);
    y1(1:neqn,1) = y2(1:neqn,1);

    t(step) = t2;
    y_0(step) = y2(1,1);
    y_1(step) = y2(2,1);

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Expected final conditions:\n' );
  fprintf ( 1, '\n' );

  if ( neqn <= 4 )
    fprintf ( 1, '  %12f', t_stop );
    for i = 1 : neqn
      fprintf ( 1, '  %12f', y_stop(i) );
    end
    fprintf ( 1, '\n' );
  elseif ( neqn <= 8 )
    fprintf ( 1, '  %12f', t_stop );
    for i = 1 : 4
      fprintf ( 1, '  %12f', y_stop(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '              ' );
    for i = 5 : neqn
      fprintf ( 1, '  %12f', y_stop(i) );
    end
    fprintf ( 1, '\n' );
  else
    y_min = min ( y_stop(1:neqn) );
    y_ave = sum ( y_stop(1:neqn) ) / neqn;
    y_max = max ( y_stop(1:neqn) );
    y_norm = sqrt ( sum ( y_stop(1:neqn).^2 ) );
    fprintf ( 1, '  %12f  %12f  %12f  %12f  %12f\n', ...
      t_stop, y_min, y_ave, y_max, y_norm );
  end

  return
end

Last modified: