I am reading "Numerical Methods for Engineers" by Chapra and Canale. In it, they've provided pseudocode for the implementation of Euler's method (for solving ordinary differential equations). Here is the pseucode:
Pseucode for implementing Euler's method
I tried implementing this code in GNU Octave, but depending on the input values, I am getting one of two errors:
error: 'ynew' undefined near line 5 column 21
error: called from
Integrator at line 5 column 9
main at line 18 column 7
I would be very grateful if you could get this program to work for me. I am actually an amateur in GNU Octave. Thank you.
Edit 1: Here is my code. For main.m:
%prompt user
y = input('Initial value of y:');
xi = input('Initial value of x:');
xf = input('Final value of x:');
dx = input('Step size:');
xout = input('Output interval:');
x = xi;
m = 0;
xpm = x;
ypm = y;
while(1)
xend = x + xout;
if xend > xf
xend = xf;
h = dx;
Integrator(x,y,h,xend);
m = m + 1;
xpm = x;
ypm = y;
if x >= xf
break;
endif
endif
end
For Integrator.m:
function Integrator(x,y,h,xend)
while(1)
if xend - x < h
h = xend - x;
Euler(x,y,h,ynew);
y = ynew;
if x >= xend
break;
endif
endif
end
endfunction
For Euler.m:
function Euler(x,y,h,ynew)
Derivs(x,y,dydx);
ynew = y + dydx * h;
x = x + h;
endfunction
For Derivs.m:
function Derivs(x,y,dydx)
dydx = -2 * x^3 + 12 * x^2 - 20 * x + 8.5;
endfunction
Edit 2: I shoud mention that the differential equation which Chapra and Canale have given as an example is:
y'(x) = -2 * x^3 + 12 * x^2 - 20 * x + 8.5
That is why the 'Derivs.m' script shows dydx to be this particular polynomial.
Your functions should look like
function [x, y] = Integrator(x,y,h,xend)
while x < xend
h = min(h, xend-x)
[x,y] = Euler(x,y,h);
end%while
end%function
as an example. Depending on what you want to do with the result, your main loop might need to collect all the results from the single steps. One variant for that is
m = 1;
xp = [x];
yp = [y];
while x < xf
[x,y] = Integrator(x,y,dx,min(xf, x+xout));
m = m+1;
xp(m) = x;
yp(m) = y;
end%while