Re: how make it faster?
- To: mathgroup at smc.vnet.net
- Subject: [mg42439] Re: [mg42429] how make it faster?
- From: Selwyn Hollis <selwynh at earthlink.net>
- Date: Tue, 8 Jul 2003 04:37:20 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Jan,
Have a look at the notebook named NumericalMethods.nb found here:
http://www.math.armstrong.edu/faculty/hollis/mmade/demos/
Basically the idea is to write a function to compute one step:
RK4step[{t_,y_}]:=
Module[{k1,k2,k3,k4}, ....; {t+h, y + h*(k1+2k2+2k3+k4)/6}]
and then use NestList:
NestList[RK4step, {t0,y0}, n]
In regard to your second question: Yes, a significant speed improvement
is possible by compiling RK4step and the function f[t,y] (where the DE
is y' = f(t,y)).
-----
Selwyn Hollis
http://www.math.armstrong.edu/faculty/hollis
On Monday, July 7, 2003, at 03:05 AM, Jan Schmedes wrote:
> Dear group,
>
> could you explain me on the following code for a simple
> Runge-Kutta-solver
> for ODE without
> stepwidth controll how it could be made faster and more
> 'mathematica-like'.
>
> RK[c_, a_, b_, func_, h_, x0_, t0_, T_] :=
> Module[{z, f, sl, s, k, i, j, ti, l, n},
> f[t_, z_] = func /. {x -> t, y -> z};
> n = Round[T/h];
> sl = Array[s, {n + 1}];
> s[1] = x0;
> ti = t0;
> For[i = 1, i <= n + 1, i++,
> j = 1;
> Do[{k[j] =
> f[ti + c[[j]]*h, s[i] + h*Sum[a[[j]][[l]]*k[l], {l, 1, j
> -
> 1}]];
> j++}, {Length[c]}];
> s[i + 1] = s[i] + h*Sum[b[[j]]*k[j], {j, 1,
> Length[c]}];
> If[i<=n-2,
> ti = ti + h;
> ,
> ti=T; (*last step*)
> ]
> ];
> Return[sl];
> ];
>
> e.g with
> (**Butcher-tableau:
> c|a
> ---
> |b
> *****)
> c = {0, 0.5};
> a = {{0, 0}, {0.5, 0}};
> b = {0, 1};
> f = (x^2 + y)/x; (*rigth side of y'=f(x,y)
> x0 = 0.;
> T = 2.;
> tau = 0.5;
> t0 = 1.;
>
> A second question is if it is possible to use Compile at some place to
> make it faster?
>
> Thank you and best regards
>
> Jan Schmedes
>
>