Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Numerical Problem

• To: mathgroup at smc.vnet.net
• Subject: [mg106982] Re: [mg106962] Numerical Problem
• From: "Tony Harker" <a.harker at ucl.ac.uk>
• Date: Sat, 30 Jan 2010 07:13:13 -0500 (EST)
• References: <201001291249.HAA29176@smc.vnet.net> <4B6318FC.3020402@wolfram.com>

]-> -----Original Message-----
]-> From: Daniel Lichtblau [mailto:danl at wolfram.com]
]-> Sent: 29 January 2010 17:21
]-> To: Tony Harker
]-> Cc: mathgroup at smc.vnet.net
]-> Subject: Re: [mg106962] Numerical Problem
]->
]-> Tony Harker wrote:
]-> > 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.7412
]-> 17},{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.7412
]-> 18},{1.3416,-0.6
]-> >
]-> 70802},{1.21415,-0.607076},{1.09881,-0.549404},{0.99442,-0.4
]-> 9721},{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.0000000
0000000000,{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.*1
0^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
]->
]-> It's not the Mathematica version that matters, it's the
]-> processor or
]-> perhaps MKL library. On one machine I get the "bad" results in both
]-> versions, and the "good" results in both versions on
]-> another amchine.

  Ah, that's interesting. I was getting different results from the different
versions on the same machine. This suggests to me that the change is in the
MKL library. But, in that all that's being done is multiplication and
addition, this seems like a rather fundamental change.

]->
]-> Here is a pared down example. First the bad result.
]->
]-> epcm[h_, y_List, f_List] := y + N[1/2* h (f.y + f.(y + h f.y))]
]->
]-> app[0] = {2, -1};
]->
]-> Do[Print[InputForm[
]->     app[j] = epcm[.1, app[j - 1], {{998, 1998}, {-999,
]-> -1999}}]]], {j,
]->    1, 6}]
]->
]-> {1.8100000000000023, -0.9049999999999955}
]-> {1.638049999944292, -0.8190249999442927}
]-> {1.4824349769833363, -0.7412173519833346}
]-> {1.3402658465670585, -0.6694638959420446}
]-> {-5.3436544706007, 5.9507302359163425}
]-> {-32138.70840490274, 32139.257808470356}
]->
]-> Here is the more desired result:
]->
]-> In[1]:= epcm[h_, y_List, f_List] := y + N[1/2* h (f.y +
]-> f.(y + h f.y))]
]->
]-> In[2]:= app[0] = {2, -1};
]->
]-> In[3]:= Do[Print[InputForm[
]->     app[j] = epcm[.1, app[j - 1], {{998, 1998}, {-999,
]-> -1999}}]]], {j,
]->    1, 6}]
]-> {1.81, -0.905}
]-> {1.63805, -0.819025}
]-> {1.48243525, -0.741217625}
]-> {1.34160390125, -0.670801950625}
]-> {1.21415153063125, -0.607075765315625}
]-> {1.0988071352212814, -0.5494035676106407}
]->
]-> In[4]:= $Version
]-> Out[4]= 8.0 for Linux x86 (32-bit) (January 28, 2010)
]->
]-> This is clearly a very ill-conditioned iteration (small changes to
]-> input, on either machine, lead to relatively larger changes
]-> in output).
]-> There may also be bad library code somewhere. To see it is ill
]-> conditioned (other than observing what is already shown),
]-> one can take
]-> derivatives with respect to the y inputs (the ones that are
]-> changing in
]-> the iteration). For simplicity we preset fh to 1/10 and f
]-> to {-999, -1999}
]->
]-> In[27]:= symbepcm[y : {y1_, y2_}] =
]->   y + 1/20*{-999, -1999}.y + {-999, -1999}.(y + 1/10*{-999,
]-> -1999}.y)
]-> Out[27]= y +
]->   1/20 {-999, -1999}.y + {-999, -1999}.(y + 1/10 {-999, -1999}.y)
]->
]-> In[33]:= D[symbepcm[{y1, y2}], y1] // N
]-> Out[33]= {298452., 298451.}
]->
]-> In[34]:= D[symbepcm[{y1, y2}], y2] // N
]-> Out[34]= {597201., 597202.}
]->
]-> These indicate that a change in input can be magnified
]-> several orders of
]-> magnitude.

  Indeed, this is to be expected, as this represents using a na=EFve numerical
scheme to solve a stiff pair of first-order differential equations in which
the required solution goes like exp(-t) and the parasitic solution is
exp(-1000 t).

]->
]-> One thing that seems to help stabilize the iteration is to
]-> use h=1/10
]-> rather than 0.1 in calling epcm[]. I doubt this is really a
]-> "fix" since,
]-> as noted, this iteration is ill conditioned.
]->

 I suppose this is picking up where in my integer-only version working in
integers is having the greatest effect.

  Tony

 
]-> Daniel Lichtblau
]-> Wolfram Research
]->
]->
]->

• References:
  • Numerical Problem
    • From: "Tony Harker" <a.harker@ucl.ac.uk>