RKF array indices error

1 次查看(过去 30 天)
Nikolaos Koutsakis
Nikolaos Koutsakis 2020-12-30
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
  1 个评论
Nikolaos Koutsakis
Nikolaos Koutsakis 2020-12-30
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.

请先登录,再进行评论。

回答(1 个)

Cris LaPierre
Cris LaPierre 2020-12-30
编辑:Cris LaPierre 2020-12-30
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)
  3 个评论
Cris LaPierre
Cris LaPierre 2020-12-30
编辑:Cris LaPierre 2020-12-30
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
Nikolaos Koutsakis 2020-12-30
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!!

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by