Trying to solve an ode using Heun's method but getting an error

Hello I am back again, I tried copying this method for Heun's method but im getting an error for some reason I cannot understand. This is annoying me.

elseif OPTION == 2

y = 1;

n = input( 'How many values should be used to calculate the numerical solution?\n')

x = linspace(0, 2, n);

h = x(2)-x(1);

for i= 1:length(x)

y_a= (x.^2+1/2+1).^2;

B(i+1) = y(i)+ h*(1+4*x(i))*sqrt(y(i));

A = x(i)+h ;

Slope_right = (1 + 4*A)*sqrt(B(1:end-1)) ;

Slope_left = (1+4*x)*sqrt(y) ;

y(i+1) = y(i)+0.5*h*Slope_left*Slope_right ;

end

plot( x, y_a,'-k', x,y(1:end-1), '-b')

legend({'Analytical','Heun'}

The error I am getting is "In an assignment A(:) = B, the number of elements in A and B must be the same."

[Jan]

This is not twitter: no # before the tags. Tahnks.

Please post the complete error message. It contains the failing line. With posting just the cause of the error, but not the location, the readers need to guess, where the problem occurs, but you have this important information on the scren already.

[Priyanshu Aryan]

this is what pops up in the command window

How many values should be used to calculate the numerical solution?

5

n =

5

In an assignment A(:) = B, the number of elements in A and B must be the same.

Answer 1

[Jan]

Slope_left = (1+4*x)*sqrt(y)

Because x is a vector, Slope_left will be a vector also. At least in the 1st iteration. In later iterations y is a vector, too, such that this line must fail.

If Slope_left is a vector, this expression is a vector also:

y(i)+0.5*h*Slope_left*Slope_right

but y(i+1) is a scalar.

[Priyanshu Aryan]

So whats the solution?

[Jan]

I do not recognize the purpose of the code. It does not look like a clean Heun algorithm. You calculate y_a, but do not use it. The intention of you "B(1:end-1)" is not clear to me also.

A bold guess is that you want to use x(i) and y(i):

Slope_left = (1+4*x(i))*sqrt(y(i));
%                  ^^^        ^^^

Search in the net for celan Heun implementations in Matlab, e.g. https://www.mathworks.com/matlabcentral/answers/349917-heun-s-method-program-code

[Priyanshu Aryan]

Let me send you the whole code maybe you can make sense of that.OPtion 1 works perfectly, however, I tried repeating the same for option 2 and 3 but didnt succedd don't know why. I havent written code for the 4th option yet.

clear, clear vars, format compact
%user inputs
OPTION = input ('Which methods do you want plotted on the graph. The options are \n 1 Euler vs. Analytical\n 2 Heun vs. Analytical \n 3 Runge-Kutta vs Analytical \n 4 All (Euler + Heun + Runge-Kutta + ode45() + Analytical \n')
if OPTION == 1 
n = input( 'How many values should be used to calculate the numerical solution?\n')
  
    y = zeros(1, n);
    y(1)= 1;
x = linspace(0, 2, n);
h = x(2)-x(1); 
for i= 1:length(x)
    y_a(i) = (x(i).^2+1/2*x(i)+1).^2;
    y(i+1) = y(i) + h*(1+4*x(i))*sqrt(y(i));
end
plot(x, y_a, '-k', x, y(1:end-1), '-b')
legend({'y_a', 'y_euler'})


elseif OPTION == 2

n = input( 'How many values should be used to calculate the numerical solution?\n')
x = linspace(0, 2, n);
y = zeros(1,n);
y(1) = 1;
h = x(2)-x(1);
B = [];
A = [];

for i= 1:length(x)
    y_a= (x.^2+1/2*x+1).^2;
    B = y(i)+ h*(1+4*x)*sqrt(y(i));
    A = x+h ;
    Slope_right = (1 + 4*A).*sqrt(B) ;
    Slope_left = (1+4*x)*sqrt(y(i)) ;
    y(i+1) = y(i)+0.5*h.*Slope_left.*Slope_right ;
    
 end
plot( x, y_a,'-k', x,y(1:end-1), '-b')
    legend({'Analytical','Heun'})





elseif OPTION == 3

n = input( 'How many values should be used to calculate the numerical solution?\n')
x = linspace(0, 2, n);
y = zeros(1,n)
    y(1) = 1;
h = x(2)-x(1);
for i= 1:length(x)
    y_a= (x.^2+1/2*x+1).^2;
k_1 = (1+4*x)*sqrt(y(i)) ; 
k_2 = (1+4*(x+0.5*h)).*sqrt(y(i)+0.5.*k_1);
k_3 = (1+4*(x+0.5*h)).*sqrt(y(i)+0.5.*k_2);
k_4 = (1+4*(x+h)).*sqrt(y(i)+ k_3);
   
y(i) = y(i)+ h*(1/6 .*(k_1)+1/3 .*(k_2)+1/3 .*(k_3)+1/6 .*(k_4));
end


    
    
    
elseif OPTION == 4
    
    y = 1;
n = input( 'How many values should be used to calculate the numerical solution?\n')
x = linspace(0, 2, n);
h = x(2)-x(1);
for i= 1:length(x)
    y_a= (x.^2+1/2+1).^2;
    y(i+1) = y(i) + h*(1+4*x(i))*sqrt(y(i));
    
    
    
    
    
    
    
    
    
    
    
    
    
   
end
end

[Jan]

The blank lines impede the reading.

This code is still not a Heuin algorithm:

for i= 1:length(x)
    y_a= (x.^2+1/2*x+1).^2;
    B = y(i)+ h*(1+4*x)*sqrt(y(i));
    A = x+h ;
    Slope_right = (1 + 4*A).*sqrt(B) ;
    Slope_left = (1+4*x)*sqrt(y(i)) ;
    y(i+1) = y(i)+0.5*h.*Slope_left.*Slope_right ;
end

1. Calculate the analyticalsolution y_a once only, not in each iteration of the loop.
2. The code would be nicer and clear, if you do not insert the function to be integrated directly in the code, but provide it as a Matlab function, e.g. "fcn".
3. Then the Heun method is easy:

fi       = fcn(x(i), y(i));
y(i + 1) = y(i) + 0.5 * h * (fi + fcn(x(i) + h, y(i) + h * fi))

[Priyanshu Aryan]

Okay thanks a lot let me try that and get back to you.