Hello! I am writing a code for a question using adaptive RK method. I encountered a problem in the while loop, where the system shows "Index Exceeds Array Bounds". I do not understand where have I gone wrong. Appreciate your help!
Note: dy_dx is a function file I wrote separately
syms dydx y xdydx = -0.6*y + 10*(exp((-(x-2)^2)/(2*(0.075^2)))) ;disp('The ODE for Q3 is:')disp(dydx)h = 0.5;e = 0.00005;% y(0) = 0.5
n = 0; i = 1;x = zeros(i+1);y = zeros(i+1);while x(i) < 4 % y1 is the est y for full step size
k1 = zeros(i); k1(i) = dy_dx(x(i),y(i)); k2 = zeros(i); k2(i) = dy_dx((x(i) + (h/2)),(y(i) + (k1(i)*h/2))); k3 = zeros(i); k3(i) = dy_dx((x(i) + (h/2)),(y(i) + (k2(i)*h/2))); k4 = zeros(i); k4(i) = dy_dx((x(i) + h),(y(i) + (k3(i)*h))); y1 = zeros(i); y1(i) = y(i) + ((h/6)*(k1(i) + 2*k2(i) + 2*k3(i) + k4(i))); % y2 is the est y for half step size
h_half = h/2; k_1 = zeros(i); k_1(i) = dy_dx(x(i),y(i)); k_2 = zeros(i); k_2(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_1(i)*h_half/2))); k_3 = zeros(i); k_3(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_2(i)*h_half/2))); k_4 = zeros(i); k_4(i) = dy_dx((x(i) + h_half),(y(i) + (k_3(i)*h_half))); y2_1 = zeros(i); y2_1(i) = y(i) + ((h_half/6)*(k_1(i) + 2*k_2(i) + 2*k_3(i) + k_4(i))); y(i) = y2_1(i); x(i) = x(i) + h_half; k_a = zeros(i); k_a(i) = dy_dx(x(i),y(i)); k_b = zeros(i); k_b(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_a(i)*h_half/2))); k_c = zeros(i); k_c(i) = dy_dx((x(i) + (h_half/2)),(y(i) + (k_b(i)*h_half/2))); k_d = zeros(i); k_d(i) = dy_dx((x(i) + h_half),(y(i) + (k_c(i)*h_half))); y2 = zeros(i); y2(i) = y(i) + ((h_half/6)*(k_a(i) + 2*k_b(i) + 2*k_c(i) + k_d(i))); del = y2(i) - y1(i); E_a = del/15; y2_new = zeros(i); y2_new(i) = y2(i) + E_a; err = abs(((y2_new(i) - y1(i))/y2_new(i))*100); if err < e break end y = y2_new(i); x = x + h_half; n = n + 1; i = i + 1; enddisp('The estimated y value is:')disp(y2_new(i-1))disp('The estimated x value is:')disp(x(i-1))disp('The number of iteration is:')disp(n)plot(x,y2_new)hold onplot(x,y);legend('Half Step Size', 'Full Step Size')xlabel('x'); ylabel('y');title('Visual Representation of the Iteration')
Best Answer