MATLAB: Adams Bashforth Multon Code not running

abmadamsbashforthmulton

Hello! I am trying to plot 3 different ABM codes with different h values as step sizes. I cannot get the script to run accurately. Any help is welcome.
Script for ABM:
function [T,Y] = abm(f,a,b,ya,h)
T = a:h:b;
Y(1) = ya;
for j = 1:min(3, length(T)-1)
    Y(j+1) = rk4(f,T(j),T(j+1),Y(j),h);
end
for j = 4:length(T)-1
    % PREDICTOR
    Y(j+1) = Y(j) + (h/24)*(55*f(T(j),Y(j)) - 59*f(T(j-1),Y(j-1)) + 37*f(T(j-2),Y(j-2)) - 9*f(T(j-3),Y(j-3)));
      % CORRECTOR
      Y(j+1) = Y(j) + (h/24)*(9*f(T(j+1),Y(j+1)) + 19*f(T(j),Y(j)) - 5*f(T(j-1),Y(j-1)) + f(T(j-2),Y(j-2)));
end
Script for RK4:
function [T,Y] = rk4(f,a,b,ya,h)
T = 2:h:10;                                         
Y = zeros(1,length(T)); 
Y(1) = ya;                                         
F_xy = @(t,y) (1/(t^2)) - ((20*y)/t);                    
for i=1:(length(T)-1)                              
    k_1 = F_xy(T(i),Y(i));
    k_2 = F_xy(T(i)+0.5*h,Y(i)+0.5*h*k_1);
    k_3 = F_xy((T(i)+0.5*h),(Y(i)+0.5*h*k_2));
    k_4 = F_xy((T(i)+h),(Y(i)+k_3*h));
      Y(i+1) = Y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  
  end
Code:
clear all
% define the problem: function f and domain
f = @(t,y) (1/(t^2)) - ((20*y)/t);
a = 2; b = 10;
% exact solution, using a fine grid
t = a:.0001:b;
y = (1./(19.*t)) - (524288./(19.*(t.^20))); % this is a vector of values, not a function
% coarse solution
h = .01;
ya = 0;
[T1,Y1]=abm(f,a,b,ya,h);
% fine solution
h = .001;
ya = 0;
[T2,Y2]=abm(f,a,b,ya,h);
% finer solution
h = .0001;
ya = 0;
[T3,Y3]=abm(f,a,b,ya,h);
plot(t,y,'k',T1,Y1,'bo-',T2,Y2,'r--',T3,Y3,'g-')
legend('Exact','h=0.01','h=0.001','h=0.0001')
title('Adams-Bashforth-Multon Method with 3 meshes')
%%ABM Relative Error
y1ex = 1./(19*T1)-524288./(19*T1.^20);
relerr1 = abs(y1ex-Y1)./(abs(y1ex)+eps);
y2ex = 1./(19*T2)-524288./(19*T2.^20);
relerr2 = abs(y2ex-Y2)./(abs(y2ex)+eps);
y3ex = 1./(19*T3)-524288./(19*T3.^20);
relerr3 = abs(y3ex-Y3)./(abs(y3ex)+eps);
plot(T1,relerr1,':',T2,relerr2,'--',T3,relerr3,'-.')

Best Answer

Look at your rk4 code. You return two variables, T and Y. The T is clearly a vector since it is assigned as
T = 2:h:10;
which is going to give a vector result unless h is negative (in which case it would be empty) or unless h exceeds +8 (in which case it would be the scalar [2]).
The rest of what you calculate in rc4 is a waste of time, as you never use the second output of rc4 in you calling code.
So you have a vector being returned from rc4, and you attempt to assign that to Y(j+1) . but j is a scalar so Y(j+1) designates a scalar location. You cannot store a vector of some 801 elements (T) inside a scalar.
You should also be paying attention to the fact that your rc4 implementation ignores its first 3 parameters. Why are you even bothering to pass them in if they are going to be ignored?
Related Question