1

I'm trying to solve the following problem using RK4:

$$y'=y^{\frac{1}{3}}, \quad y(0)=0, \quad t\in[0,6].$$

The exact solution is $y(t)=\sqrt{\left( \frac{2t}{3} \right)^{3}}$. I wrote the following code:

int main(){

  double a = 0;
  double b = 6;

  double t_0 = 0;
  double y_0 = 0; 

  int n = 10000;
  double h = (b-a)/n;

  double t[n+1];
  double y_RK[n+1];

  t[0]=t_0;
  y_RK[0]=y_0;

  double k1,k2,k3,k4;

  for(int i=0; i<n; ++i){

    k1 = f(t[i],y_RK[i]);
    k2 = f(t[i] + h/2, y_RK[i] + (h/2)*k1);
    k3 = f(t[i] + h/2, y_RK[i] + (h/2)*k2);
    k4 = f(t[i] + h,y_RK[i]+h*k3);

    t[i+1]=t[i]+h;
    y_RK[i+1] = y_RK[i] + (h/6)*(k1 + 2*(k2+k3) + k4);

  }

}

However, the result I'm getting is this: Red = RK4 result, Green = actual solution

Why isn't it working? Sorry for remaking the question, I realised I was doing it wrong previously, but now I rewrote the code and it still doesn't converge. I'm not sure if this exercise is meant to illustrate that RK4 doesn't always work (the question says "investigate"), or if I'm making a mistake?

Spine Feast
  • 4,770
  • I found another wonderful tidbit. A friend told me this problem must have a mistake in it, since it's unable to "take off" from zero. In my case it did, because in C++, apparently $0^{\frac{1}{3}} = 1$. – Spine Feast Oct 26 '14 at 20:20

1 Answers1

4

There is no hope it can work when starting at $y = 0$ (and any $t$), because the Cauchy-Lipschitz conditions are not satisfied on the $x$-axis. Indeed, the function $y\mapsto y^{1/3}$ is not locally lipschitzian at $y = 0$, which is visible since it has an infinite derivative at this point.

Another consequence is that the equation has actually an infinity of solutions passing through the point $(y=0,t=0)$. Indeed, consider the function which is $0$ from $0$ to some real $a$, and then a solution similar to yours, but starting at $(y=0,t=a)$. This makes a perfectly acceptable solution of the equation.

See http://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem for more informations.

user187552
  • 56