Question: Runge-Kutta of order 2 does not evaluate

Hi all,

Suppose that you want to solve x'(t)=f(x,t)   x(t0)=x0

The theory show that:

x(t+h)=x(t)+1/2*h*(F1+F2) where F1=f(x,t) and F2=f(x+F1,t+h)

SO I made procedure:

> RK2 := proc (fct, xinit, tinit, a, b, N)
 local h, x, t, F1, F2, i;
 x[0] := xinit; t[0] := tinit; h := (b-a)/N;
 F1 := unapply(fct, [x, t]);
 F2 := unapply(subs({x = x+F1(x, t), t = t+h}, F1(x, t)), [x, t]);
 for i from 0 to N do
 x[i+1] := evalf(x[i]+(1/2)*h*(F1(x[i], t[i])+F2(x[i], t[i])));
 t[i+1] := evalf(t[i]+(i+1)*h)
 end do; 
[seq([t[i], x[i]], i = 0 .. N)]
 end proc;

If we take, for the sake of an exemple,  f(x,t)=1-t+x/t

fct:=1-t+x/t and x(1)=0

Then

> P := RK2(fct, 0, 1, 1, 2, 10);

You will see that for a strange reason,  there is some "x" and "t" that are not evaluate.

But it suppose to work since if I do it manually, it work!

> fct := 1-t-x/t;

> F1 := unapply(fct, [x, t]);

> F2 := unapply(subs({t = t+h, x = x+F1(x, t)}, F1(x, t)), [x, t]);

> h := 1/10; x[0] := 0; t[0] := 1;

> x[1] := evalf(x[0]+(1/2)*h*(F1(x[0], t[0])+F2(x[0], t[0])));

                           x[1] := -0.005000000000
> t[1] := 1+h;

                                         11
                                 t[1] := --
                                         10
> x[2] := evalf(x[1]+(1/2)*h*(F1(x[1], t[1])+F2(x[1], t[1])));

                           x[2] := -0.01558712122
and so on...

What am I missing out?

Thanks in advance.

Mario
 

Please Wait...