Numerical Problem
- To: mathgroup at smc.vnet.net
- Subject: [mg106962] Numerical Problem
- From: "Tony Harker" <a.harker at ucl.ac.uk>
- Date: Fri, 29 Jan 2010 07:49:12 -0500 (EST)
I expect I'm missing something in the documentation, but I can't see why
this simple example of Euler's predictor-corrector applied to a stiff pair
of equations should behave differently in Version 7 and Version 5.
========================================================================
$Version
epcm[h_,y_List,f_List]:=y+0.5 h (f.y+f.(y+h f.y))
fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
app1=fepcm[.1,0.,1,{2.,-1.},{{998,1998},{-999,-1999}}]
7.0 for Microsoft Windows (32-bit) (February 18, 2009)
{{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48243,-0.741217},{1.34027,-0.
669464},{-5.34365,5.95073},{-32138.7,32139.3},{-1.57517*10^8,1.57517*10^8},{
-7.71992*10^11,7.71992*10^11},{-3.78353*10^15,3.78353*10^15},{-1.85431*10^19
,1.85431*10^19}}
============================================================================
$Version
epcm[h_,y_List,f_List]:=y+0.5 h (f.y+f.(y+h f.y))
fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
app1=fepcm[.1,0.,1,{2.,-1.},{{998,1998},{-999,-1999}}]
5.2 for Microsoft Windows (June 20, 2005)
{{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48244,-0.741218},{
1.3416,-0.670802},{1.21415,-0.607076},{1.09881,-0.549404},{0.99442,-0.\
49721},{0.899951,-0.449975},{0.814455,-0.407228},{0.737082,-0.368541}}
============================================================================
It's clear that the version 5 result is more acceptable, as we can do
everything in integers and convert to reals at the end:
$Version
epcm[h_,y_List,f_List]:=y+ h (f.y+f.(y+h f.y))/2
fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
N[app1=fepcm[1/10,0,1,{2,-1},{{998,1998},{-999,-1999}}]]
7.0 for Microsoft Windows (32-bit) (February 18, 2009)
{{2.,-1.},{1.81,-0.905},{1.63805,-0.819025},{1.48244,-0.741218},{1.3416,-0.6
70802},{1.21415,-0.607076},{1.09881,-0.549404},{0.99442,-0.49721},{0.899951,
-0.449975},{0.814455,-0.407228},{0.737082,-0.368541}}
============================================================================
What I don't understand is why this discrepancy between Versions 5 and 7
occurs, and whether I can get version 7 to behave in the same way as 5.
For something even more exciting, there's
============================================================================
$Version
epcm[h_,y_List,f_List]:=y+0.500000000000000000 h (f.y+f.(y+h f.y))
fepcm[h_,t0_,tmax_,y0_List,f_List]:=NestList[epcm[h,#,f]&,y0,10]
app1=fepcm[.100000000000000000,0.00000000000000000,1.00000000000000000,{2.00
000000000000000,-1.00000000000000000},{{998.000000000000000,1998.00000000000
0000},{-999.00000000000000000,-1999.00000000000000000}}]
7.0 for Microsoft Windows (32-bit) (February 18, 2009)
{{2.0000000000000000,-1.0000000000000000},{1.810000000000,-0.905000000000},{
1.6380500,-0.8190250},{1.482,-0.741},{0.*10^1,0.*10^1},{0.*10^6,0.*10^6},{0.
*10^10,0.*10^10},{0.*10^15,0.*10^15},{0.*10^20,0.*10^20},{0.*10^24,0.*10^24}
,{0.*10^29,0.*10^29}}
============================================================================
Tony Harker