Implementing Euler's Method in GNU Octave

Question

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:

The program doesn't give any output at all. I have to press 'Ctrl + C' in order to break execution.
The program gives this message:

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.

Without your code, we can not see how it compares to the given code and where the error might come from. And if you interpreted the interface of the pseudo code correctly. — Lutz Lehmann, Apr 07 '20 at 11:26
That is better. Next adapt the code from a Pascal/C like calling convention to a matlab/octave/scilab like function call structure, like `function [xo,yo] = Integrator(x,y,h,xnew)` and then in the main loop use `[x,y] = Integrator(x,y,dx,min(xend, x+dx))` or similar. Also, look up on how to create arrays for `xp,yp` and how to fill them with values. // For your question, add the parameter that you tested that give the 2 reported results. // My scilab installation does not work with the Java version, so at the moment no complete error lists. — Lutz Lehmann, Apr 07 '20 at 12:49
I can no longer tell where I got the idea from that this is a scilab question... — Lutz Lehmann, Apr 07 '20 at 16:22

score 2 · Answer 1 · answered Apr 07 '20 at 15:38

Here is my final code. It has four different M-files:

main.m

%prompt the user
y = input('Initial value of y:');
x = input('Initial value of x:');
xf = input('Final value of x:');
dx = input('Step size:');
xout = dx;

%boring calculations
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

%plot the final result
plot(xp,yp);
title('Solution using Euler Method');
ylabel('Dependent variable (y)');
xlabel('Independent variable (x)');
grid on;

Integrator.m

%This function takes in 4 inputs (x,y,h,xend) and returns 2 outputs [x,y]
function [x,y] = Integrator(x,y,h,xend)
  while x < xend
    h = min(h, xend-x);
    [x,y] = Euler(x,y,h);
  end
endfunction

Euler.m

%This function takes in 3 inputs (x,y,h) and returns 2 outputs [x,ynew] 
function [x,ynew] = Euler(x,y,h)
   dydx = Derivs(x,y);
   ynew = y + dydx * h;
   x = x + h;
endfunction

Derivs.m

%This function takes in 2 inputs (x,y) and returns 1 output [dydx]
function [dydx] = Derivs(x,y)
    dydx = -2 * x^3 + 12 * x^2 - 20 * x + 8.5;
endfunction

Yes, that is almost identical to the other parts I did not post. Note that you can put subroutines also in the same .m file as the main function of some method (some decoration required). Where that makes sense depends on the size of the procedures, and how re-usable the method is. Here one could sensibly combine 2. and 3. if in the function parameters also the handle for the derivatives function is passed. — Lutz Lehmann, Apr 07 '20 at 16:20

Lutz Lehmann · Accepted Answer · 2020-04-07T16:23:31.387

1

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

edited Apr 07 '20 at 16:23

answered Apr 07 '20 at 13:19

Lutz Lehmann

25,219
2
22
51

Thank you so much! Your solution works perfectly, and I also got to learn about how functions work in Octave (I'm a complete amateur). BTW, do you know any good books from which to learn about "basic numerical methods" for GNU Octave? I'm an electronics engineering student, so my programming skills are not too strong. – Anuj Manoj Shah Apr 07 '20 at 15:29
remember that Octave is built to be Matlab code compatible. So for resources you can look at what's been don for both programs over the years. MIT Opencourseware has a numerical methods course with a Matlab/Octave primer - https://ocw.mit.edu/courses/nuclear-engineering/22-15-essential-numerical-methods-fall-2014/index.htm – Nick J Apr 08 '20 at 13:42
here are the 3-part tutorial videos: https://ocw.mit.edu/courses/nuclear-engineering/22-15-essential-numerical-methods-fall-2014/tutorial-videos/ – Nick J Apr 08 '20 at 13:42

Implementing Euler's Method in GNU Octave

2 Answers2