Adams Bashforth Multon Code not running

Question

Kaylene Widdoes 2016-2-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/269827-adams-bashforth-multon-code-not-running

编辑： Meysam Mahooti 2018-12-30

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,'-.')

2 个评论
显示无隐藏无

Jan 2016-2-24

Please explain why you assume that you cannot get the code to run "accurately".

Kaylene Widdoes 2016-2-24

It won't run at all

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2016-2-24

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/269827-adams-bashforth-multon-code-not-running#answer_211104

在 MATLAB Online 中打开

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?

2 个评论
显示无隐藏无

Kaylene Widdoes 2016-2-24

So how should I fix the problem while still being able to use different values of h I designate at multiple points?

Walter Roberson 2016-2-25

You would start by documenting your code. What are the inputs to each routine? What are the outputs from each routine? What algorithm is being used to do the calculations? If the inputs are of such-and-such a size, what size will the outputs be?

Then you would proceed through your code and see if your code matches the documentation. Do you use all of the inputs? If not then what point is there in passing in that input? Do you assign values to all of the outputs? In the places that you call the routines, are all of the inputs passed? In the places that you call the routines, is a variable assigned for each of the outputs? Has the calling code assigned the right amount of storage for the expected output, not too big and not too small?

Having a clear understanding of what the inputs are and what the outputs are, and ensuring that the inputs and outputs match those expectations, goes a long way towards getting any routine working.

请先登录，再进行评论。