RE: multistep iterative methods
- To: mathgroup at smc.vnet.net
- Subject: [mg40342] RE: [mg40324] multistep iterative methods
- From: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
- Date: Wed, 2 Apr 2003 04:35:40 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
>-----Original Message----- >From: Selwyn Hollis [mailto:hollisse at mail.armstrong.edu] To: mathgroup at smc.vnet.net >Sent: Tuesday, April 01, 2003 11:52 AM >To: mathgroup at smc.vnet.net >Subject: [mg40342] [mg40324] multistep iterative methods > > >I'd like to throw this out as a challenge to the group: What's >the most >efficient way to implement in Mathematica an explicit multistep >iterative method such as, say, the 4-step Adams-Bashforth method for >solving y' = f(t,y): > >y[k+1]:= y[k] + (h/24)*(55*f[k] - 59*f[k-1] + 37*f[k-2] - 9*f[k-3]) > >where y[0], y[1], y[2], y[3] are "given," and f[i] denotes f[t0 +i*h, >y[i]]. The desired output would be the list > >{y[0], y[1], y[2], ... , y[n]}. > >A suitable toy problem is > >y' = -2t*y^2, y(0) = 1, > >with h = 0.01, n = 1000 (?), and the starting values taken from the >exact solution y = 1/(1+t^2): > >y[0]=1, y[1] = 0.9999, y[2] = 0.9996, y[3] = .999101. > >Thanks in advance. > >------- >Selwyn Hollis > > > > Selwyn, I think, I'd better learn that from you! Well, here my (naive) attempt: In[2]:= f[t_, y_] := -2 t y^2 In[12]:= t0 = 0.; h = 0.01; y[0] = 1; y[1] = 0.9999; y[2] = 0.9996; y[3] = .999101; In[13]:= dd = {-9, 37, -59, 55}*h/24 In[14]:= initstep := (fk = f[t0 + h*#, y[#]] & /@ Range[0, 3]; t = t0 + 3 h) In[17]:= step = (ys = # + dd.fk; fk = RotateLeft[ReplacePart[fk, f[t += h, ys], 1]]; ys) & In[21]:= yex[t_] := 1/(1 + t^2) ------------- In[24]:= (ylist = Join[{y[0], y[1], y[2]}, initstep; NestList[step, y[3], 997]]); // Timing Out[24]= {0.15 Second, Null} ------------- In[25]:= xylist = Transpose[{h*Range[0, 1000], ylist}]; In[26]:= gapp = ListPlot[xylist, PlotJoined -> True, PlotStyle -> {Thickness[.003], Hue[.5, 1, .7]}] In[27]:= gex = Plot[yex[t], {t, 0, 1000*h}, DisplayFunction -> Identity] In[28]:= Show[gapp, gex, PlotRange -> All] -- Hartmut Wolf