#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

`help/text/ITER`:=TEXT(
`HELP FOR: numerical solution of 1st order d.e. for any iterative scheme`,
` `,
`CALLING SEQUENCE`,
`   ITER(args)`,
` `,
`SYNOPSIS`,
`-input required to be of the form`,
`    ITER(de,ischemex,ischemey,nsteps,stepsiz,xc1,yc1,xc2,yc2,..)`,
` where de is the differential equation,`,
`       ischemex is the iterative scheme for x (similarly for y)`,
` used to evaluate y(x) at nsteps intervals of size stepsiz`,
` from the starting conditions xc1,yc1,xc2,yc2`,
`-two starting conditions are required for the leapfrog method`,
` compared to only one for other iterative schemes`,
`-gfun and ischeme are functions e.g.`,
`       gfun:=(x,y)->x^2;`,
`       ischemey:=(x,y,stepsiz)->y+stepsiz*gfun(x,y);`,
`       ischemex:=(x,stepsiz)->x+stepsiz;`,
`-stepsiz is included as`,
` a variable parameter as a half interval is used in the mid-point`,
` method for example`,
`-the actual solution is shown for comparison`,
` `,
`EXAMPLE:`,
`Given the differential equation`,
`      dy/dx=x^2`,
`and the starting condition (x=1,y=1) obtain a numerical approximation`,
`for y when x=7 using the Euler method`,
` `,
`To proceed define the parameters to be used in IRATE:`,
`      de:=diff(y(x),x)=x^2;`,
`      gfun:=(x,y)->x^2;`,
`      ischemey:=(x,y,stepsiz)->y+stepsiz*gfun(x,y);`,
`      ischemex:=(x,stepsiz)->x+stepsiz;`,
`Now call the procedure:`,
`      ITER(de,ischemex,ischemey,6,1,1,1);`,
`The 4th and 5th parameters can be changed to improve the approximation`,
` `,
`To use a better iterative scheme simply redefine ischemex and`,
`ischemey e.g. for the leapfrog method:`,
`      ischemey:=(y1,x2,y2,stepsiz)->y1+2*stepsiz*gfun(x2,y2);`,
`      ischemex:=(x,stepsiz)->x+2*stepsiz;`,
`The procedure call then requires an additional starting condition`,
`      ITER(de,ischemex,ischemey,6,1,1,1,2,2);`
):

#............................................................................
#note of iterative relations that can be used
#1.Euler
#   y[n+1]=y[n]+hf(x[n],y[n])
#2.Mid-point
#   y[n+1]=y[n]+hf(x[n]+h/2,y[n]+h/2f(x[n],y[n]))
#3.Heun
#   y[n+1]=y[n]+h/2(f(x[n],y[n])+f(x[n+1],y[n]+hf(x[n],y[n]))
#4.Leapfrog
#   y[n+2]=y[n]+2hf(x[n+1],y[n+1])
#____________________________________________________________________________

ITER:=proc()
#parameters should be entered in this order
#      de,ischemex,ischemey,nstep,stepsiz,xc1,yc1,xc2,yc2,vecout,vecin,plotop
#where the last parameters are multiple starting conditions
#for a scheme like LEAPFROG
#nargs=9 for Euler and similar, nargs=11 for LEAPFROG
       local xc,yc,l,L,i,j,p,q,r,ans,xi,yi,
             ischemex,ischemey,nstp,choix1,choix2,choix3,choix4;
       global vecout, vecin, g,const,xs,xb,gfun,grad;
#vecout stores the output coords yc in a global vector
#vecin provides input for a de as discrete values of a function of
#the dependent variable, particularly for integrating the output
#from the stochastic de/random walk

#parameters must be passed in a certain order, but search all 
#parameters to discover which options are required i.e leapfrog
#or something to do with the stochastic problem

#after setting choices make sure parameters are correctly assigned to
#variables. 1st 7 variables are the default - Euler type method with
#no options
nstp:=op(2,args[4])+1;
xc:=array(1..nstp);yc:=array(1..nstp);

ischemex:=args[2];ischemey:=args[3];

#starting conditions
xc[1]:=op(2,args[6]);yc[1]:=op(2,args[7]);

for i from 1 to nargs do
  if op(1,args[i])=`y2` then print(`Ah Leapfrog!`);
                             choix1:=1;
                             xc[2]:=xc[1]+op(2,args[5]);#print(xc[2]);
                             yc[2]:=op(2,args[i]); fi;
  if op(1,args[i])=`noplot` then print(`no analytical soln/plot`);
                                 choix2:=1 fi; 
  if op(1,args[i])=`numeric` then print(`dsolve(  ,numeric)`);
                                 choix2:=2 fi;
  if op(1,args[i])=`dataout` then print(`output stored in vector`);
                                  choix3:=1; 
                                  vecout:=op(2,args[i]) fi;
  if op(1,args[i])=`datain` then print(`input data from vector`);
                                 choix4:=1; 
                                 vecin:=op(2,args[i]) fi;
od;



if choix3=1 then vecout[1]:=yc[1]; fi;


#calculate the right limit to this loop for the no of steps args[4]
if choix1<>1 then
#Euler or similar
  for i from 1 to op(2,args[4]) do
#   print(i);
    xc[i+1]:=ischemex(xc[i],op(2,args[5]));#print(");
    yc[i+1]:=ischemey(xc[i],yc[i],op(2,args[5]));#print(");
    if choix4=1 then yc[i+1]:=subs(op(2,args[1])=vecin[i],yc[i+1]); fi;
    if choix3=1 then vecout[i+1]:=yc[i+1]; fi;
  od;
else
#leapfrog
  for i from 2 to op(2,args[4]) do
    xc[i+1]:=ischemex(xc[i-1],op(2,args[5]));#print(xc[i]);
    yc[i+1]:=ischemey(yc[i-1],xc[i],yc[i],op(2,args[5]));
  od;
#print(xc);
fi;

L:=[];
#loop to construct plot list
for j from 1 to i-1 do
  l:=[[xc[j],0],[xc[j],Re(yc[j])],[xc[j+1],Re(yc[j+1])],[xc[j+1],0]];
  L:=[op(L),op(l)]; #print(");
od;

if choix2=1 then
  plot(L,x=xc[1]..xc[i],style=LINE,title=`numerical solution`);
elif choix2=2 then
#solve de numerically and plot
  xi:=op(2,args[6]);yi:=op(2,args[7]);
  ans:=dsolve({args[1],y(xi)=yi},y(x),numeric);
  q:=odeplot(ans,[x,y(x)],xc[1]..xc[i],numpoints=25);
  p:=plot(L,x=xc[1]..xc[i],style=LINE,title=`numerical solution`);

  display({p,q}); 
else
#solve de and plot
  xi:=op(2,args[6]);yi:=op(2,args[7]);
  ans:=dsolve({args[1],y(xi)=yi},y(x));
  q:=Re(op(2,ans));
  r:=plot(q,x=xc[1]..xc[i]);
  p:=plot(L,x=xc[1]..xc[i],style=LINE,title=`numerical solution`);

  display({p,r}); 
fi;

end:


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

`help/text/QUOT`:=TEXT(
`HELP FOR:quotients procedure called QUOT`,
` `,
`CALLING SEQUENCE`,
`   QUOT(args)`,
` `,
`SYNOPSIS`,
`-input required to be of the form`,
`   QUOT(f,x0,h,x1)`,
` where f is the function whose gradient is to be estimated`,
` at the point x0 using an interval of length h starting at x1`,
`-the approximate gradient is compared with the actual gradient`,
` at x0. `,
` `,
`EXAMPLE:`,
`To obtain an estimate of the gradient of a function y`,
`    y:=cos(x);`,
`at x=1 using a foward difference quotient with interval length 0.5`,
`use QUOT:`,
`    QUOT(y,1,0.5,1);`,
`If a backward difference is required the last parameter defining the`,
`start of the interval should be 0.5 and for a centred difference it`,
`should be 0.75.`
):

#____________________________________________________________________________

QUOT:=proc(y,x0in,hin,x1in)
      local h,x0,x1,xl,xh,x2,f,q,Q,l,xb,xs;
      global  grad,g,const,gfun, gf;  
x0:=op(2,x0in);h:=op(2,hin);x1:=op(2,x1in);
#calculate the difference quotient (foward, backward or centred)
#depending upon the input parameters x1 and x2
xl:=x0-3*h;
xh:=x0+3*h;
q:=plot(y,x=xl..xh,title=`difference quotients`);
Q:={q};

x2:=x1+h;
f:=makeproc(y,x);
#grad:=(f(x2)-f(x1))/h;
l:=[[x1,0],[x1,f(x1)],[x2,f(x2)],[x2,0]];
q:=plot(l,style=LINE);
Q:=Q union {q};

l:=[[x0,0],[x0,f(x0)]];
q:=plot(l,style=LINE);
Q:=Q union {q};

#show the actual gradient aswell
g:=makeproc(diff(y,x),x);
grad:=g(x0);
const:=f(x0)-grad*x0;
gfun:=x->grad*x+const;
xs:=x0-2*h;xb:=x0+2*h;
l:=[[xs,gfun(xs)],[xb,gfun(xb)]];
q:=plot(l,style=LINE);
Q:=Q union {q};

display(Q);

end:


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#with(plots):#with(student):