在GNU Octave中实施Euler方法

Question

我正在阅读Chapra和Canale撰写的“工程师的数值方法”。在其中，他们为实现Euler方法（用于求解常微分方程）提供了伪代码。这是伪代码：

Pseucode for implementing Euler's method

我曾尝试在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的业余爱好者。谢谢。

Edit 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

Edit 2：我应该提到Chapra和Canale给出的微分方程为：

y'（x）= -2 * x ^ 3 + 12 * x ^ 2-20 * x + 8.5

这就是'Derivs.m'脚本显示dydx是该特定多项式的原因。

Answer 1

您的函数应该看起来像

  function [x, y] = Integrator(x,y,h,xend)
    while x < xend
      h = min(h, xend-x)
      [x,y] = Euler(x,y,h);
    end
  endfunction

例如。根据您要对结果执行的操作，您的主循环可能需要从单个步骤中收集所有结果。一种方法是

  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

Answer 2

这是我的最终代码。它具有四个不同的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;

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

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

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