What do I need to change to get my Improved Euler's to match my Euler's?

My Euler's code is displaying the correct plot but the Improved Euler's shows the line as increasing infinitely. It should be closer to the Euler's plot. it is the same for my Runge-Kutta equation

% Euler's Part 2
clear x;
clear y;

LD = 4; % Last digit of student ID
SD = 1; % second-to-last digit of student ID
A = sqrt(10+LD);
B = 0.05*(1+LD+SD);
V = B % diameter of valve

V = 0.3000

y0 = 2*A % y(0)

y0 = 7.4833

g = 9.8

g = 9.8000

xi = 0;
h = 0.5;
xf = 650;

n = (xf - xi)/h;
x(1) = xi;
y(1) = y0;

% -((V^2)/(4))*sqrt((2*g)/(y^3))

for j=1:n
    x(j+1)=h+x(j);
    y(j+1)=y(j)+h*(-((V^2)/(4))*sqrt((2*g)/((y(j))^3)));
end

figure(4)
plot(x,y)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

xlabel('t')
ylabel('y(t)')
title('Euler')
hold on
grid on
% Improved Euler's Part 2
clear x
clear y

LD = 4; % Last digit of student ID
SD = 1; % second-to-last digit of student ID
A = sqrt(10+LD);
B = 0.05*(1+LD+SD);
V = B % diameter of valve

V = 0.3000

y0 = 2*A % y(0)

y0 = 7.4833

g = 9.8

g = 9.8000

xi = 0;
h = 0.5;
xf = 650;

n = (xf - xi)/h;
x(1) = xi;
y(1) = y0;

% -((V^2)/(4))*sqrt((2*g)/(y^3))

for j=1:n
    x(j+1)=h+x(j);
    f(j)=-((V^2)/(4))*sqrt((2*g)/((y(j))^3));
    y(j+1)=y(j) + h*f(j)/2 + (h/2)*(h -((V^2)/(4))*sqrt((2*g)/((y(j))^3))-(h*f(j)));
end
figure(4)
plot(x,y)
xlabel('t')
ylabel('y(t)')
title('Improved Euler')
hold on
grid on


% RK4 Part 2
clear x
clear y

LD = 4; % Last digit of student ID
SD = 1; % second-to-last digit of student ID
A = sqrt(10+LD);
B = 0.05*(1+LD+SD);
V = B % diameter of valve

V = 0.3000

y0 = 2*A % y(0)

y0 = 7.4833

g = 9.8

g = 9.8000

xi = 0;
h = 0.5;
xf = 650;

n = (xf - xi)/h;
x(1) = xi;
y(1) = y0;

% -((V^2)/(4))*sqrt((2*g)/(y^3))

for j=1:n
    x(j+1)=h+x(j);
    f(j)= -((V^2)/(4))*sqrt((2*g)/((y(j))^3));

    k1(j)=f(j);
    k2(j)= h/2 -((V^2)/(4))*sqrt((2*g)/((y(j))^3)) -(h*k1(j)/2);
    k3(j)= h/2 -((V^2)/(4))*sqrt((2*g)/((y(j))^3)) -(h*k2(j)/2);
    k4(j)= h -((V^2)/(4))*sqrt((2*g)/((y(j))^3)) -(h*k3(j));

    y(j+1)= y(j)+(h*k1(j)/6)+(h*k2(j)/3)+(h*k3(j)/3)+(h*k4(j)/6);
end

figure(4)
plot(x,y)
xlabel('t')
ylabel('y(t)')
title('Runge-Kutta')

hold on
grid on

Answer 1

[James Tursa]

Your current code is hard to read and debug because you have multiple variations of this expression scattered throughout your code:

(-((V^2)/(4))*sqrt((2*g)/((y(j))^3)))

Rather than trying to track down your errors with your current code, I suggest you rewrite everything using a function handle for the derivative. E.g.,

f = @(x,y) -V^2/4*sqrt(2*g/y^3)

Then in your subsequent lines of code call this function handle as f(x(j),y(j)) etc. instead of repeating the expression over and over again. This will greatly facilitate examining your code for errors, and also it allows you to call the MATLAB supplied functions such as ode45( ) for direct comparison, as well as comparing your code to pseudo-code that can be found on Wikipedia et al.

Try to do this, then repost everything for us to look over. I can see what look like obvious errors in your calculation and usage of the k's, for instance, but rather than work with what you currently have it would be best for a rewrite first. E.g., your Euler integration line becomes simply:

y(j+1) = y(j) + h * f(x(j),y(j));

Your Improved Euler would involve the equivalent of

h * ( f(x(j),y(j)) + f(x(j+1),y(j+1)) ) / 2 

on the rhs, as long as you calculate the Euler step y(j+1) first (maybe store it in a temporary variable). These types of expressions in your code will be much easier to debug for errors than what you currently have.

The RK4 code becomes much simpler and would look something like this:

k1 = f(x(j)      , y(j)         );
k2 = f(x(j) + h/2, y(j) + k1*h/2);
k3 = f(x(j) + h/2, y(j) + k2*h/2);
etc.

[Jessica]

I get this error:

% Euler's Part 2
clear x;
clear y;

LD = 4; % Last digit of student ID
SD = 1; % second-to-last digit of student ID
A = sqrt(10+LD);
B = 0.05*(1+LD+SD);
V = B % diameter of valve

V = 0.3000

y0 = 2*A % y(0)

y0 = 7.4833

g = 9.8

g = 9.8000

xi = 0;
h = 0.5;
xf = 650;

n = (xf - xi)/h;
x(1) = xi;
y(1) = y0;

% -((V^2)/(4))*sqrt((2*g)/(y^3))
f = @(x,y) -V^2/4*sqrt(2*g/y^3)

f = function_handle with value:
    @(x,y)-V^2/4*sqrt(2*g/y^3)

for j=1:n
    x(j+1)=h+x(j);
    f(j)=f(x(j), y(j));
    y(j+1)=y(j) + h * ( f(x(j),y(j)) + f(x(j+1),y(j+1)) ) / 2 ;
end

Conversion to function_handle from double is not possible.

figure(4)
plot(x,y)
xlabel('t')
ylabel('y(t)')
title('Euler')
hold on
grid on

[James Tursa]

This line needs to be an Euler step as I mentioned before:

f(j)=f(x(j), y(j))

So change it to the same code as your Euler step that I showed you:

y(j+1) = y(j) + h * f(x(j),y(j));

Then this y(j+1) gets used in the next line to calculate the Improved Euler step.

The reason your current code is getting an error is because you have the function handle f(etc) on both the lhs and rhs. It should only be on the rhs.