在 GNU Octave 中实现欧拉方法

Question

我正在阅读 Chapra 和 Canale 的 "Numerical Methods for Engineers"。在其中，他们提供了用于实施欧拉方法（用于求解常微分方程）的伪代码。这是伪代码：

Pseucode for implementing Euler's method

我尝试在 GNU Octave 中实现此代码，但根据输入值，我收到以下两个错误之一：

该程序根本不提供任何输出。我必须按 'Ctrl + C' 才能中断执行。
程序给出了这条信息：

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：我应该提到 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%while
  end%function

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

  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

Answer 2

这是我的最终代码。它有四个不同的 M 文件：

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

在 GNU Octave 中实现欧拉方法

Implementing Euler's Method in GNU Octave

octave

numerical-methods

ode

differential-equations

numerical-analysis