Runge-Kutta method, cannot produce the correct graphs

The code should output, three graphs evaluating the 2nd order Runge-Kutta method for varying step size. However, currently the output is a graph with three horizontal lines, that should be curving upwards and tending towards the actual solution. I am confused as to what i have done wrong, any advice would be great.

Answer 1

[James Tursa]

You need to index x in your calculations. E.g.,

        q1 = f(x(i),y(i));
        q2 = f(x(i)+h, y(i)+h*q1);

You should initialize y(1), not y(2):

        y(1) = 1;

And I would set up your indexing and for-loop this way:

    x = (2:h:4);
    no_points = numel(x);
    for i = 1:no_points-1  

And plot the entire x and y vectors

        plot(x,y,'o-')

[Ben Bawtree]

thank you! that makes sense, unfortunately the intial data we have been given is y(2)=1, so i have to stick to that, does this affect the output?

[James Tursa]

You are misunderstanding the meaning of "y(2) = 1" in your instructions. That simply means that when x=2 the corresponding value of y is 1. It does not mean that the index value used for your y vector in your code has to be 2. The index values you use are different from the actual x and y values. In the code I have suggested you have this setup

x(1) = 2

y(1) = 1

That is, when the initial value of x is 2, the corresponding initial value of y is 1. That is what you want. The (1) indexing is just that ... only indexing. It is not the values. Bottom line is that the y(1) = 1 statement is correct and matches your instructions because it is paired with the x(1)=2 value.

[Ben Bawtree]

ah thank you! i appreciate you explaining, cheers

Answer 2

[William Rose]

@Ben Bawtree,

See attached code. You assigned a value to y(2) initially:

y(2)=1;

This creates a vector of length 2, where y(1)=0 by default, and y(2)=1. This is kind of strange because i t means y(2) goes with x(1), etc. So I changed that by asisgning a value to y(1):

y(1)=1;

Then I plotted the x vector versus the y vector after doing the 2nd order RK integration. I changed the outer loop variable from h=[.5 .2 .1] to j=1:3. This allowed me to have a value j that I could use to specify which third of the figure to use for each plot. As a result, I had to define h as a 3 element vector before entering the outer loop. Then, inside the loop, I refer to h(j) wherever I need the value of h.

See code and see plot of results.

%BenBawtreeCode.m    WCRose 20210311
clear all; close all; clc;
f = @(x,y) 5*y-3;
%next line creates a vector of length 2 and assigns a value to the second element
%y(2) = 1;           
%I think you actualy want to assign a value to first element of y
y(1)=1;
h=[.5 .2 .1];
figure;
for j = 1:3 
    % define the number of points to be evaluated;
    no_points = 2/h(j)+1;
    x = (2:h(j):4);
%    for i = 2:no_points  
    for i = 1:no_points-1  
        % RK2 method 
        q1 = f(x,y(i));
        q2 = f(x+h(j), y(i)+h(j)*q1);
        y(i+1) = y(i) + h(j)*0.5*(q1+q2);
    end
    disp(['h = ',num2str(h(j)), ', y(4)=', num2str(y(i+1))]);
    subplot(3,1,j);             %plot in a specific sub-third of the figure
    plot(x,y,'ro-');            %plot all the x values versus all the y values
    grid on;
    ylabel('Y');        
    titlestr=sprintf('h=%.1f',h(j));
    title(titlestr);            %plot titleh
    if j==3, xlabel('X'); end   %X axis label on bottom plot
end

[William Rose]

@Ben Bawtree,

Notice that your function

f(x,y)=5*y-3

defines the derivative. Therefore the diffential equation you are solving is

dy/dx = 5y-3.

This is easily solved with standard calculus. The solution is

y(x) = K*exp(5*x)+0.6.

The unknown constant K is determined by applying your knowledge of th initial conditions. IN this case, we know from the code that y=1 when x=2. We use this informaiton to solve for K:

K=0.4/exp(10).

You can add the analytic solution to your plot. You will see that when h gets small (h<=.001), the RK2 solution approaches the analytic solution, which, to repeat, is

y(x) = K*exp(5*x)+0.6.

Answer 3

[William Rose]

@Ben Bawtree,

This code plots the solutions with different values of h on a single plot, and it adds the analytic solution.

%BenBawtreeCode.m    WCRose 20210311
clear;
f = @(x,y) 5*y-3;       %dy/dx
%Analytic solution to  dy/dx=5y-3 is K*exp(5x)+0.6.
%Since y=1 at x=2, it must be that K=0.4/exp(10).

y(1)=1;
h=[.5 .2 .1 1e-3];
K=0.4/exp(10);
figure; 
hold on;
plotstylestr=["ro-","go-","bo-","mo-","k.-"];   %We will use this below
for j = 1:length(h) 
    x = (2:h(j):4);
    for i = 1:length(x)-1  
        % RK2 method 
        q1 = f(x,y(i));
        q2 = f(x+h(j), y(i)+h(j)*q1);
        y(i+1) = y(i) + h(j)*0.5*(q1+q2);
    end
    disp(['h = ',num2str(h(j)), ', y(x=4)=', num2str(y(i+1))]);
    plot(x,y,plotstylestr(j));    %plot all the x values versus all the y values
end 
plot(x,K*exp(5*x)+.6,plotstylestr(end));
grid on;
ylabel('Y');  xlabel('X');  %labels on axes
%Make sure next line matches values in h(). 
legend('h=0.5','h=0.2','h=0.1','h=0.001','Analytic');

Plot generated by code above:

[Ben Bawtree]

Thank you for your answer! Very helpful, will have a look over my code