Question

I have to write procedure for euler's method to differential equations in Pascal. I have this:

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.

Code of procedure Euler() is ok - I copied it from book. When I compile it on Lazarus I received error:

Project project1.exe raised exception class 'External: SIGSEGV'

So I trying run my app directly from folder via windows console - and then I got:

An unhandled exception occurred at $00000000 : EAccessViolation : Access violation $00000000 $00401A78 main, line 123 of project1.lpr

I trying compile it in another compiler (Dev-Pascal) - it's working but shows nothing. I don't understand what happen. Anyone have any idea what's wrong?

Was it helpful?

Solution

Correct your code

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);

Variable n_f was not be initialized by function.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top