1

我必须为欧拉方法编写帕斯卡微分方程的程序。我有这个：

program Euler_proc;
  {$mode objfpc}{$H+}
uses
  {$IFDEF UNIX}{$IFDEF UseCThreads}
  cthreads,
  {$ENDIF}{$ENDIF}
  Classes
  { you can add units after this };

{$IFDEF WINDOWS}{$R project1.rc}{$ENDIF}

{THIS PART I JUST COPIED FROM BOOK: }

type fxy = function (x,y : Extended) : Extended;
function f(x,y:Extended): Extended; far;
begin
f:= x+x*y+y+1;
end;

procedure Euler (var x0 : Extended;
                 x1     : Extended;
                 var y  : Extended;
                 f      : fxy;
                 mh,eps : Extended;
                 var h  : Extended;
                 var fl : Boolean;
                 var st : Integer);

const c=3.333333333333333333e-1;
var   h1,h2,g,g1,x,xh,y1,y2,yh : Extended;
      hend,tend                : Boolean;

begin
  st:=0;
  if fl
    then begin
           h:=x1-x0;
           fl:=false
         end;
  if h<=0
    then st:=1
    else begin
           tend:=true;
           hend:=true;
           x:=x0;
           repeat
             g:=h/2;
             g1:=g/2;
             yh:=y+g*f(x,y);
             xh:=x+g;
             y1:=y+h*f(xh,yh);
             yh:=y+g1*f(x,y);
             xh:=x+g1;
             y2:=y+g*f(xh,yh);
             xh:=x+g;
             yh:=y2+g1*f(xh,y2);
             y2:=y2+g*f(xh,yh);
             xh:=ln((1+c)*abs(y1-y2)/eps);
             h1:=h/exp(c*xh);
             if h1<=c*h
               then if 2*h1<mh
                      then begin
                             st:=2;
                             hend:=false;
                             y:=y2;
                             x:=x+h;
                             h2:=h
                           end
                      else h:=2*h1
               else if h1>=h
                      then if x+2*h<=x1
                             then h:=2*h
                             else begin
                                    h2:=h;
                                    h:=x1-x;
                                    x:=x+h;
                                    tend:=false
                                  end
                      else begin
                             y:=y2;
                             x:=x+h;
                             if x+2*h1<x1
                               then h:=2*h1
                               else begin
                                      h:=x1-x;
                                      x:=x+h;
                                      h2:=2*h1;
                                      if h2>x1-x0
                                        then h2:=x1-x0;
                                      tend:=false
                                    end
                           end
           until not tend or not hend;
           if hend
             then begin
                    g:=h/2;
                    xh:=x+g;
                    yh:=y+g*f(x,y);
                    y:=y+h*f(xh,yh)
                  end;
           x0:=x;
           h:=h2
         end
end;
   {THIS IS MY CODE:}

   var n_st:Integer;
       n_h :Extended;
       n_f :fxy;
       n_y :Extended;
       n_x0:Extended;
       n_fl:Boolean;
       n_x1:Extended;
       n_mh:Extended;
       n_esp:Extended;

BEGIN
//Euler(x0,x1,y,f,mh,esp,h,fl,st);
n_x0:=-1;
n_y:=1;
n_fl:=True;
n_x1:=1;
n_mh:=1e-4;
n_esp:=1e-8;

Euler(n_x0,n_x1,n_y,n_f,n_mh,n_esp,n_h,n_fl,n_st);

writeln(n_x0);
writeln(n_y);
writeln(n_h);
writeln(n_fl);

readln();
readln();
END.

程序代码 Euler() 没问题 - 我从书中复制了它。当我在 Lazarus 上编译它时，我收到了错误：

项目 project1.exe 引发异常类“外部：SIGSEGV”

所以我尝试通过 Windows 控制台直接从文件夹运行我的应用程序 - 然后我得到：

在 $00000000 处发生未处理的异常：EAccessViolation：访问冲突 $00000000 $00401A78 main，project1.lpr 的第 123 行

我尝试在另一个编译器（Dev-Pascal）中编译它 - 它正在工作但什么也没显示。我不明白发生了什么。有人知道出了什么问题吗？

4

1 回答 1

3

更正您的代码

n_esp:=1e-8;

n_f := @f; // Add this line

Euler(n_x0,n_x1,n_y,n_f,n_mh,n_esp,n_h,n_fl,n_st);

writeln(n_x0);

变量n_f未由函数初始化。

于 2013-06-01T09:48:16.150 回答