2

I'm trying to determine how many steps it takes for the Runge-Kutta Method ("RK4") to be within 0.01% of the exact solution of an ordinary differential equation. I'm comparing this to the Euler method. Both should result in a straight line on a loglog plot. My Euler solution seems to be correct but I am getting a curved solution for RK. They are based on the same code so I am completely confused about the problem.

EDIT: Sorry for removing the wikipedia links. It wouldn't let me keep more than one link since I'm a new user. However, both methods are detailed on Wikipedia similarly to my implementation.

Should someone wish to tackle my problem, the code is below and graphs are in this Word file on dropbox.com. Yes this is a homework problem; I'm posting this because I'd like to understand what is wrong in my thought process.

f = @(x,y) x+y; %this is the eqn (the part after the @(t,y)

This is my RK4 code:

k1=@(x,y) h*f(x,y);
k2=@(x,y) h*f(x+1/2*h,y+1/2*k1(x,y));
k3=@(x,y) h*f(x+1/2*h,y+1/2*k2(x,y));
k4=@(x,y) h*f(x+h,y+k3(x,y));

clear y x exact i
x(1)=0;
y(1)=2;
xn=1;
x0=0;

exact=3*exp(xn)-xn-1; %exact solution at x=1
%# Evaluate RK4 with 1 step for x=0...1
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
RK4(N)=y(i+1);  %# store result of RK4 in vector RK4(# of steps)
E_RK4(N)=-(RK4(N)-exact)/exact*100;%keep track of %error for each N
Nsteps_RK4(N)=N;


%# repeat for increasing N until error is less than 0.01%
while -(RK4(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
        x(i+1)=x(i)+h;
    end
    RK4(N)=y(i+1);
    Nsteps_RK4(N)=N;
    E_RK4(N)=-(RK4(N)-exact)/exact*100; %# keep track of %error for each N
end

Nsteps_RK4(end);

This is my Euler code: 

%# Evaluate Euler with 1 step for x=0...1
clear y x i
x(1)=0;
y(1)=2;
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)= y(i)+h*f(x(i),y(i));
Euler(N)=y(i+1); %# store result of Euler in vector Euler(# of steps)
E_Euler(N)=-(Euler(N)-exact)/exact*100;%# keep track of %error for each N
Nsteps_Euler(N)=N;
%# repeat for increasing N until error is less than 0.01%
while -(Euler(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)= y(i)+h*f(x(i),y(i));   
        x(i+1)=x(i)+h;
    end
    Euler(N)=y(i+1);
    Nsteps_Euler(N)=N;
    E_Euler(N)=-(Euler(N)-exact)/exact*100; %# keep track of %error for each N
end
eykanal
  • 26,437
  • 19
  • 82
  • 113
Dominik
  • 782
  • 7
  • 27
  • 1
    How is the function `f` defined? – stardt Oct 09 '11 at 22:34
  • Hi, as stardt said you should tell us which `f`you are supposed to use. Also, can you make sure the Wikipedia article links I added correctly reflect the version of RK4 you are using? – Jonas Heidelberg Oct 09 '11 at 22:44
  • There are many problems with the code you posted. `h` is used in the anonymous functions before it is defined. You should use `abs()` in the `while` conditions because I don't see how you can guarantee that the estimate will always be less than the exact solution. You access the wrong element of the `y` vector after computing the values in the `for` loop. The line after `end` should be `RK4(N) = y(N+1)`. The same mistakes are in the euler code. – stardt Oct 09 '11 at 22:50
  • Also, changing `h` in the `while` loop does not affect the `h` used in the k1...k4 functions. Consider:`>> h = 1; >> g = @(x) x+h; >> g(0) ans = 1 >> h = 2; >> g(0) ans = 1` – stardt Oct 09 '11 at 22:58
  • I'm wrong about the indexing and for loop variable scope issue. You still have a problem with changing `h`. The code as it is posted here will not run because `h` is undefined. To see this, use `clear all` instead of clearing only specific variables. The `h` in your `k` functions is whatever value `h` had when you defined those functions. I suggest using N=N*2 for the explicit Euler code to get the results much faster. – stardt Oct 10 '11 at 01:04
  • N=i at the end of each for loop so both produced the same results. I did also change the while loop to use abs() as it's cleaner. Using h in my function was the main problem. I'm not sure I understand though, why the functions of k1...k4 would not use the new h from each iteration. I'll try to rewrite the code with k1...k4 not as functions but as scalars changing with each loop. – Dominik Oct 10 '11 at 01:35
  • with a different `f` you would need to change the `while` logic because the approximate answer could be greater than exact answer. I don't see how avoiding `abs` makes the code "cleaner". With `abs` the same code works with any `f`. – stardt Oct 10 '11 at 01:59

1 Answers1

1

See: http://www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html#f4-71621

Variables specified in the body of the expression [...] must have a value assigned to them at the time you construct an anonymous function that uses them. Upon construction, MATLAB captures the current value for each variable specified in the body of that function. The function will continue to associate this value with the variable even if the value should change in the workspace or go out of scope.

The value of h in k1...k4 remains constant even though you change h in the while loop.

One solution is to add h to the anonymous functions:

k1=@(x,y,h) h*f(x,y);
stardt
  • 1,179
  • 1
  • 9
  • 14
  • Thank you! That explains how those functions work. And adding h as a variable of the function is much easier than rewriting the lops with a scalar h. I could use one clarification, I ran the code with both `k1=@(x,y,h) h*f(x,y);` and `k1=@(x,y,h) h*f(x(i),y(i));` and both worked. MATLAB seemed to recognize i as an index instead of a scalar (since I hadn't defined i yet). Is this a correct interpretation or did I miss a difference in the output? – Dominik Oct 10 '11 at 06:08
  • @Dominik My guess is that you had `i` being 1 when you defined the anonymous function. If `i` were undefined you would have obtained an error when calling `k1()`. I like to use `clear all` at the top of my scripts to make sure I don't have unexpected variables. – stardt Oct 10 '11 at 15:17