2

I need to use Mathematica to solve the following DE by the fourth-order Runge-Kutta method and find $y(10)$: $$y'(x)=\frac{1}{2y}, y(0)=1$$ The following is the code I have so far (with step size h=0.5):

Clear[x, y, n, h, S1, S2, S3, S4, i];
y[0] = 1; y[1] = Sqrt[2]; h = 0.5; n = 20;
Do[x[i] = 0 + (i - 1)*h, {i, 1, n + 1}]

y'[x_] = f[x_, y_] = 1/(2*y);
Do[{S1 = f[x[i], y[i]],
  S2 = f[x[i] + h/2, y[i] + h*S1/2],
  S3 = f[x[i] + h/2, y[i] + h*S2/2],
  S4 = f[x[i] + h, y[i] + h*S3],
  y[i + 1] = y[i] + (S1 + 2 S2 + 2 S3 + S4)*h/6}, {i, 1, n}]
Do[Print[x[i], " ", y[i]], {i, 1, n + 1}]

This code results in an approximation of $y(10)=3.4641$, which is fine; however, when I change the step size and the corresponding $n$ value I still get $y(10)=3.4641$ regardless of the value of the step size, which doesn't make sense to me. I'm almost certain that the smaller the step size one chooses, the better the approximation should be, so we shouldn't observe the same value in error. I'm not sure how to proceed from here, but I'm sure there is an error in my code somewhere. Any help is greatly appreciated.

linearalgebrathrowaway
  • 123
  • 2

1 Answers1

1

You are suppose to give the initial condition y[1]=1 not y[0]=1 because your loop starts from i=1.

Your this condition y[1]=Sqrt[2] doesn't make any sense.

If you start your loop from i=0 and ends at i=n then you can use y[0]=1?

For h=0.01, I get y[10]=3.31632. Comparing this with the results given by NDSolve,

Clear[x, y];
sol = NDSolve[{y'[x] == 1/(2*y[x]), y[0] == 1}, y[x], {x, 0, 10}];
y[x] /. sol /. x -> 10

3.31662

Clear[x, y, S1, S2];
a = 0; b = 10;
y[0] = 1; h = 0.01; n = (b - a)/h;
Do[x[i] = 0.0 + (i)*h, {i, 0, n}]

f[x_, y_] = 1/(2*y);
Do[{S1 = f[x[i], y[i]],
  S2 = f[x[i] + h, y[i] + h*S1],
  S3 = f[x[i] + h/2, y[i] + h*S2/2],
  S4 = f[x[i] + h, y[i] + h*S3],
  y[i + 1] = y[i] + (S1 + 2 S2 + 2 S3 + S4)*h/6}, {i, 0, n}]
Do[Print[x[i], "       ", y[i]], {i, 0, n}];
zhk
  • 11,939
  • 1
  • 22
  • 38