RKF array indices error

Hello and Merry Christmas to all,

i am trying to solve an ode using RKF with 6th order accuracy.

The form of the ode is y'(t) = f(t,y(t)) and specifically is :

y' = sqrt(y), with initial values y(0) = 1 and the exact solution is y(t) = 1/4(t + 2)^2

i am using this code that i get from mathworks documentation about RKF method

clc; close all; clear all;

%% Inputs and Declaration
h = 0.1; % solution stepsize
x = 0:h:1; % input vector

% Function declaration
f = @(x, y) sqrt(y);

% Initial conditions
y_RKF(0) = 1;
gamma = [16/135, 0, 6656/12825, 28561/56430, -9/50, 2/55]; % vector related to RKF procedure

%% Core: Runge Kutta 6th order procedure
for i = 1:length(x)-1
k1 = h*f(x(i), y_RKF(i));                                                       k1 = double(k1);
k2 = h*f(x(i)+h/4, y_RKF(i)+k1/4);                                              k2 = double(k2);
k3 = h*f(x(i)+3/8*h, y_RKF(i)+3/32*k1+9/32*k2);                                 k3 = double(k3);
k4 = h*f(x(i)+12/13*h, y_RKF(i)+1932/2197*k1-7200/2197*k2+7296/2197*k3);        k4 = double(k4);
k5 = h*f(x(i)+h, y_RKF(i)+439/216*k1-8*k2+3680/513*k3-845/4104*k4);             k5 = double(k5);
k6 = h*f(x(i)+h/2, y_RKF(i)-8/27*k1+2*k2-3544/2565*k3+1859/4104*k4-11/40*k5);   k6 = double(k6);
K = [k1, k2, k3, k4, k5, k6];
y_RKF(i+1) = y_RKF(i) + sum(K.*gamma); % new solution
end

%% Analytical Exact Solution
y_exact = 1/4 .* (x + 2) .^ 2;
error_RKF = y_exact - y_RKF; % error between 

But i get this error :

Array indices must be positive integers or logical values.

Error in rkf4_5 (line 13)

y_RKF(0) = 1;

Any help why this happens? Thanks in advance

[Nikolaos Koutsakis]

The y_RKF(0) = 1 indicates the initial conditions of the ode or something entirely else? cause if i change it to y_RKF(1) = 1 it works but it doesn't output the right solution cause the initial conditions are wrong. Thanks in advance.

Answer 1

[Cris LaPierre]

y_RKF(0) = 1 creates an error because MATLAB does not use 0-based indexing. All indexing must use positive integers (1, 2, 3...) or logical indexing.

What your initial condition y(0)=1 means is y=1 when t=0. This corresponds to and index of 1 because x(1)=0.

Last thought - your error is currently the elementwise difference between corresponding values in y_exact and y_RKF. You probably want the total error, so sum the result.

error_RKF = sum(y_exact - y_RKF)

[Nikolaos Koutsakis]

Thanks so much for your quick reply! But, i think that i didn't express my question properly. So please let me rephrase it one more time. The initial conditions for this ode is y(0)=1, so how can i implement it to this script? Since the

y_RKF(0) = 1

doesn't work. I change it to

y_RKF(1) = 1

like the example of the documentation for the RKF method, it works properly with no errors, but obviously it's not right cause the initial conditions are wrong. Any ideas how it works? Or how can i implement the initial condition y(0)=1 in the same way like the example? Here is the example https://www.mathworks.com/matlabcentral/fileexchange/62581-runge-kutta-fehlberg-rkf-numerical-method.

Also, thanks for the advice about the sum of the two solutions but i need the difference in each step of the method, so the sum doesn't really help. But still, i am really really gratefull!

[Cris LaPierre]

You are confusing function input and indexing. Here, y is a function of t. In your code, you define t as the variable x.

h = 0.1; % solution stepsize
x = 0:h:1; % input vector

So when you say you want y(0)=1, you mean the value of y at t=0 to be 1, or

x and y are variables, each containing values. You want to set the initial condition, or the value of y when t=0, to 1. The first value of x is 0, so the corresponding value in y (the first value) should be set to 1.

x(1)

ans = 0

% this works
y(1)=1

y = 1

% This does not
y(0)=1

Array indices must be positive integers or logical values.

[Nikolaos Koutsakis]

Thanks so much!! I get now what you mean and you made me realize how it works! Thanks you so much again for your time!!

Merry Christmas and Happy New Year!!