Реализация метода Эйлера в GNU Octave

Question

Я читаю «Численные методы для инженеров» Чапры и Канале. В нем они предоставили псевдокод для реализации метода Эйлера (для решения обыкновенных дифференциальных уравнений). Вот псевдокод:

Псевдок для реализации метода Эйлера

Я пытался реализовать этот код в GNU Octave, но в зависимости от входных значений я получаю одно из двух ошибки:

1. Программа вообще не выдает никаких результатов. Я должен нажать 'Ctrl + C', чтобы прервать выполнение.
2. Программа выдает следующее сообщение:

error: 'ynew' undefined near line 5 column 21
error: called from
    Integrator at line 5 column 9
    main at line 18 column 7

Я был бы очень признателен, если бы вы могли получить эта программа для меня. Я на самом деле любитель в GNU Octave. Спасибо.

Редактировать 1 : Вот мой код. Для 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

Для 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

Для Euler.m:

function Euler(x,y,h,ynew)
   Derivs(x,y,dydx);
   ynew = y + dydx * h;
   x = x + h;
endfunction

Для Derivs.m:

function Derivs(x,y,dydx)
    dydx = -2 * x^3 + 12 * x^2 - 20 * x + 8.5;
endfunction

Редактировать 2 : Я должен упомянуть, что дифференциальное уравнение, которое Чапра и Канале привели в качестве примера:

y '(x) = -2 * x ^ 3 + 12 * x ^ 2 - 20 * x + 8,5

Именно поэтому скрипт Derivs.m показывает, что dydx является именно этим полиномом.

Answer 1

Вот мой окончательный код. Он имеет четыре разных M-файла:

1. 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

Answer 2

Ваши функции должны выглядеть как

  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

в качестве примера. В зависимости от того, что вы хотите сделать с результатом, вашему основному l oop может потребоваться собрать все результаты за один шаг. Один вариант для этого -

  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