CPS615 Homework 7

- Ho-Yen Chang



Assignment 7, Due date: November 6, 1995

For this problem set, write HPF programs to investigate the solution of ordinary differential equations. Turn in program listings and output Programs should have sufficient comments to allow readers to understand them with ease. A discussion of the results from your output should also be included.

1a. Use the Runge-Kutta (Order Four) algorithm to solve the initial value problem:

     y' = y - (t** 2) + 1
     0 <= t ,= 2
     y(0) = 0.5

Solutions should be generated for the time points t= 0.0, 0.2, 0.4...2.0


Solution :

Program rk1.hpf

The result output rk1.out

******************************************************************************

1b. Provide empirical evidence that the Runge-Kutta (order four) algorithm exhibits fourth-order convergence. Use the following strategy. Solve the above initial value problem using the 0.001, 0.005, 0.010, 0.020 and 0.040 time steps. Assume that the results for the 0.001 time step are correct. Calculate the cumulative errors for each of the 0.005, 0.010, 0.020 and 0.040 time step series relative to the 0.001 time step. Finally,compare the 0.005, 0.010, 0.020 and 0.040 time step cumulative errors to determine the convergence order.


Solution :

Program rk2.hpf

The result output rk2.out

The complete result reference rk2.ref

 #######################################################
 #                                                     #
 # Using Runge -Kutta (Order Four) algorithm to solve  #
 #    y' = y - (t**2) + 1                              #
 #    0 <= t <=2 , y(0) = 0.5                          #
 #    time steps h = 0.001, 0.005, 0.010, 0.020, 0.040 #
 #                                                     #
 # Assume that the result for 0.001 time step are      #
 # correct and calculate the cummulative errors for    #
 # each of the time steps series relative to the 0.001 #
 #                                                     #
 #######################################################
  
 ************ when h=0.001, N=2000 ************
 This is to be used as a standard one 
 the cumulative errors relative to itself h=0.001 is 0.0
  
  
 ************** when h=0.005, N=400 **************
  
  
 the cumulative errors relative to h = 0.001 is  7.0634681437198310E-009
  
 ************** when h=0.01, N=200 **************
  
  
 the cumulative errors relative to h = 0.001 is  5.6670162784655531E-008
  
 ************** when h=0.02, N=100 **************
  
  
 the cumulative errors relative to h = 0.001 is  4.5491147537823906E-007
  
 ************** when h=0.04, N=50 **************
  
  
 the cumulative errors relative to h = 0.001 is  3.6631778929097081E-006

"cerror" = the cumulative errors relative to h = 0.001 If we make every cerror/(h**i), i= 1, 2, 3, 4, 5 , and then we will get the following results :


 ------- the result of that cerror/(h**1) ---------
 c2[1] (h=0.005) ====> 1.4126936287439662E-06
 c3[1] (h=0.01)  ====> 5.6670162784655531E-06
 c4[1] (h=0.02)  ====> 2.2745573768911953E-05
 c5[1] (h=0.04)  ====> 9.1579447322742702E-05

 ------- the result of cerror/(h**2) ---------
 c2[2] (h=0.005) ====> 2.8253872574879324E-04
 c3[2] (h=0.01)  ====> 5.6670162784655531E-04
 c4[2] (h=0.02)  ====> 1.1372786884455977E-03
 c5[2] (h=0.04)  ====> 2.2894861830685675E-03

 ------- the result of cerror/(h**3) ---------
 c2[3] (h=0.005) ====> 5.6507745149758641E-02
 c3[3] (h=0.01)  ====> 5.6670162784655524E-02
 c4[3] (h=0.02)  ====> 5.6863934422279876E-02
 c5[3] (h=0.04)  ====> 5.7237154576714182E-02

 ------- the result of cerror/(h**4) ---------
 c2[4] (h=0.005) ====> 1.1301549029951728E+01
 c3[4] (h=0.01)  ====> 5.6670162784655522E+00
 c4[4] (h=0.02)  ====> 2.8431967211139937E+00
 c5[4] (h=0.04)  ====> 1.4309288644178546E+00

 ------- the result of cerror/(h**5) ---------
 c2[5] (h=0.005) ====> 2.2603098059903455E+03
 c3[5] (h=0.01)  ====> 5.6670162784655520E+02
 c4[5] (h=0.02)  ====> 1.4215983605569969E+02
 c5[5] (h=0.04)  ====> 3.5773221610446363E+01

By the observation of the results above, when i=3, the value of cerror/(h**3) in all the time steps seems to converge to a same constant, that is, the cumulative errors are very likely to converge when order=3. Thus, we can guess the convergence order equal to 3. However, it may depends on other factors, because we just use the h=0.001 as our standard value of comparision, and it does not mean that this value is the exact correct value.

******************************************************************************

2a .Use the Adasm Fourth-Order Predictor-Corrector algorithm to solve the initial value problem:

     y' = y - (t** 2) + 1
     0 <= t <= 0.5

Solutions should be generated for the time points t = 0.0, 0.2, 0.4 ...2.0


Solution :

Program ad1.hpf

The result output ad1.out

******************************************************************************

2b. Provide empirical evidence that the Adams Fourth-Order Predictor-Corrector algorithm exhibits fourth-order convergence. Use the following strategy. Solve the above initial value problem using the 0.001, 0.005, 0.010, 0.020 and 0.040 time step series relative to the 0.001 time step. Finally, compare the 0.005, 0.010, 0.020 and 0.040 time step cumulative errors to determine the convergence order.


Solution :

Program ad2.hpf

The result output ad2.out

The complete result reference ad2.ref

 #######################################################
 #                                                     #
 # Using Adams Fourth-Order Predictor-Corrector        #
 # algorithm to solve                                  #
 #    y' = y - (t**2) + 1                              #
 #    0 <= t <=2 , y(0) = 0.5                          #
 #    time steps h = 0.001, 0.005, 0.010, 0.020, 0.040 #
 #                                                     #
 # Assume that the result for 0.001 time step are      #
 # correct and calculate the cummulative errors for    #
 # each of the time steps series relative to the 0.001 #
 #                                                     #
 #######################################################
  
 ************ when h=0.001, N=2000 ************
 This is to be used as a standard one 
 the cumulative errors relative to itself h=0.001 is 0.0
  
 ************** when h=0.005, N=400 **************
  
 the cumulative errors relative to h = 0.001 is  1.3499270368200200E-008
  
 ************** when h=0.01, N=200 **************
  
 the cumulative errors relative to h = 0.001 is  1.0568627406559727E-007
  
 ************** when h=0.02, N=100 **************
  
 the cumulative errors relative to h = 0.001 is  8.0781187472567240E-007
  
 ************** when h=0.04, N=50 **************
  
 the cumulative errors relative to h = 0.001 is  5.9083192336384371E-006

"cerror" = the cumulative errors relative to h = 0.001 If we make every cerror/(h**i), i= 1, 2, 3, 4, 5 , and then we will get the following results :


 ------- the result of cerror/(h**1) ---------
 c2[1] (h=0.005) ====> 2.6998540736400400E-06
 c3[1] (h=0.01)  ====> 1.0568627406559727E-05
 c4[1] (h=0.02)  ====> 4.0390593736283620E-05
 c5[1] (h=0.04)  ====> 1.4770798084096093E-04

 ------- the result of cerror/(h**2) ---------
 c2[2] (h=0.005) ====> 5.3997081472800801E-04
 c3[2] (h=0.01)  ====> 1.0568627406559727E-03
 c4[2] (h=0.02)  ====> 2.0195296868141810E-03
 c5[2] (h=0.04)  ====> 3.6926995210240232E-03

 ------- the result of cerror/(h**3) ---------
 c2[3] (h=0.005) ====> 1.0799416294560160E-01
 c3[3] (h=0.01)  ====> 1.0568627406559727E-01
 c4[3] (h=0.02)  ====> 1.0097648434070905E-01
 c5[3] (h=0.04)  ====> 9.2317488025600580E-02

 ------- the result of cerror/(h**4) ---------
 c2[4] (h=0.005) ====> 2.1598832589120317E+01
 c3[4] (h=0.01)  ====> 1.0568627406559727E+01
 c4[4] (h=0.02)  ====> 5.0488242170354525E+00
 c5[4] (h=0.04)  ====> 2.3079372006400143E+00

 ------- the result of cerror/(h**5) ---------
 c2[5] (h=0.005) ====> 4.3197665178240641E+03
 c3[5] (h=0.01)  ====> 1.0568627406559726E+03
 c4[5] (h=0.02)  ====> 2.5244121085177261E+02
 c5[5] (h=0.04)  ====> 5.7698430016000358E+01

By the observation of the results above, when i=3, the value of cerror/(h**3) in all the time steps seems to converge to a same constant, that is, the cumulative errors are very likely to converge when order=3. Thus, we can guess the convergence order equal to 3. However, as the discussion of 1b., it may depends on other factors, because we just use the h=0.001 as our standard value of comparision, and it does not mean that this value is the exact correct value.

******************************************************************************


References: Burden, R.L. and Faires, J.D. (1993). Numerical Analysis, Fifth Edition, Prindle, Weber and Schmidt, Boston, pp.259-261, 275-280.


Ho-Yen Chang
B17-8 Slocum Heights, Syracuse, NY 13210
Tel (315) 442-7239
hchang03@top.cis.syr.edu