# File ohjelmat.mpl  spring 1999
# This text file contains Maple codes developed, used and/or found 
# to be used/useful/usable in mat-1.174, Computational PDE course in
# Spring 1999
# Download this file into your directory and read into your Maple session
# by the Maple-command:
# > read("path/ohjelmat.mpl");
#
interface(showassumed=2);  # Better place:  $HOME/.mapleinit 
# 
#    Green's function for Lu=-u''
 G:=(x,y)->piecewise(y>=0 and y <=x,y*(1-x),x<=y and y<=1,x*(1-y));
#### Green's function for Sturm Liouville.
 Grf:=(x,s)->piecewise(s>=0 and s <=x,v1(x)*v0(s),x<=s and
 s<=1,v0(x)*v1(s));
#########################
# Example:
# f:=x->x;
# assume(x>0);additionally(x<1);
# u:=int(G(x,y)*f(y),y=0..1);
#
## Code taken from Betounes, there seem to be some problems with it.
# 
# DYNAMIC DEplots in MAPLE
#  
# To view the heat in the plate flowing dynamically from the top edge of
# the plate, where it is hottest, toward the respective boundaries,
# which are held at zero degrees, we have to augment MAPLE's graphing
# capabilities with the following code. This creates a procedure called
# DDEplot  which functions much like MAPLE's  DEplot procedure (that you
# loaded with the with(DEtools) command above). You should execute the
# following code now in order to get this procedure defined.
# 
 DDEplot:=proc(deqs,vars,trange,initset,xrange,yrange,stsize)local
 de,ini,M,tr,xr,yr,T0,T1,H,a,b,c,d,N,ps,j,x0,y0,f,f1,f2,n,pL,k;
 de:=deqs;ini:=initset; M:=nops(ini);tr:=trange;xr:=xrange;yr:=yrange;
 T0:=op(1,op(2,tr));T1:=op(2,op(2,tr));H:=op(2,stsize);
 a:=op(1,op(2,xr));b:=op(2,op(2,xr));c:=op(1,op(2,yr));d:=op(2,op(2,yr)
 );
 N:=eval(floor((T1-T0)/H)); 
 ps:=array(1..N,1..M);
 for j to M do x0:=op(2,op(j,ini));y0:=op(3,op(j,ini));
 f:=dsolve({op(1,de),op(2,de),x(T0)=x0,y(T0)=y0},{x(t),y(t)},type=numeric,output=listprocedure,method=classical);
 f1:=subs(f,x(t));f2:=subs(f,y(t));
 for n to N do ps[n,j]:=seq([f1(k*H),f2(k*H)],k=1..n);od;od;
 pL:=eval(seq([seq(CURVES([ps[n,j]],COLOUR(RGB,0,0,0)),j=1..M)],n=1..N)
 );
 PLOT(ANIMATE(pL),AXESLABELS(x,y),VIEW(a..b,c..d)); end:
 
 
# The procedure DDEplot has the same syntax as DEplot  but it recognizes
# fewer options. This means that you can take any DEplot command , add a
# D to it, and then execute it. The extra options (if any) will be
# ignored by DDEplot. More importantly, the output from DDEplot will be
# an animation, which you can manipulate with all the controls in the
# animation window. Thus for example:
# 
# CD EXERCISES (7b)
# (2) Execute the following to get an animation of part of the static
# picture in Figure 7.12. CAUTION: this can take several minutes.
# DDEplot([deqs],[x,y],t=0..2,{init2},x=0..1,y=0..1,stepsize=.05);
# 
# Within the animation window, you can slow the animation down, move it
# ahead one frame at a time, or reverse the direction of flow. It is
# also interesting to choose the points style for the animation. In the
# points style, note that the distance between the points gets smaller
# as time goes on. Does this seem reasonable based on the mountain slope
# analogy. Estimate the average speeds at various parts on the ski
# trails (see the Note below):
# 
# NOTE: In DDEplot, the number of frames in the animation is the same as
# the number of time steps used in the numerical algorithm for the
# static DEplot. Thus if the time interval is   t=t[0]...t[1] and the
# stepsize is  h, then the number of time steps (frames) is  N =
# (t[1]-t[0])/h  (or the greatest integer less than this number).
# 
# 
#################
# Maple-oppaan ohjelmat
#Tässä on listaus ohjelmakoodeista, jotka voi hakea alla olevasta 
#URL:stä. Koneella olevaan tiedostoon tulee lisäyksiä ja parannuksia
#nopemmin kuin kirjan sivulle. 

## Tämä elää ja on jo kehittynyt kirjan versioista.
##########################################################
## Viimeinen päivitys 18.3. 1999 ##
##########################################################
#Jos haluat saada nämä funktiot käyttöösi, niin hae ne www-selaimesi
#{\tt File}-valikon {\tt Save as} valinnan kautta.
#Kaikki proseduurit saat käyttöön \maple-istunnossa komentamalla
#
#> read(`ohjelmat.mpl`);

#Alla tiedoston kommenteissa kerrotaan, miten polku annetaan \win-ympäristössä,
#\unix:ssa se tehdään, kuten aina.


#Tämänhetkinen ohjelmalistaus on puutteellinen ja keskeneräinen.
#Tiedostoa ylläpidetään ja sen voi ladata www-sivun kautta.
#
# http://www.math.hut.fi/~apiola/maple/opas/ohjelmat.mpl
#=========================================================
#
# Download this file into your working directory. 
# In your Maple session do
# > read("path/ohjelmat.mpl");
# In pc the path syntax is of the form
# > read("c:\\windows\\mymaple\\ohjelmat.mpl");
#
# Yleiskäyttöisiä apufunktioita
#  
# Vector-list-conversions:
#############################
v2l:=vec->convert(vec,list):
l2v:=lis->convert(lis,vector):

# Liitetään kaksi vektoria peräkkäin:
vappend:=(v1,v2)->vector([op(convert(v1,list)),op(convert(v2,list))]):

m2l:=M->map(op,convert(M,listlist)):  # rowwise "ravel" the 
                                      # matrix into a long list                # Compare to Matlab's M(:) which does it columnwise
#
#### Diskretointi vrt. Matlab:n linspace #############
# 
linspace:=proc() local i,n,a,b; 
          a:=args[1];b:=args[2];
          if nargs=2 then n:=10 else n:=args[3] fi;
          [seq(a+i*(b-a)/(n-1),i=0..n-1) ] end:
# Esim: x10:=linspace(0,1); x20:=linspace(0,1,20);
#
# Kts. 5.3 Diskretointi s. 117. Vaihdoimme oletusarvon 10:een, Matlab:n 100;n
# sijasta.
#
#### ITERAATIOT ################
#
 porras := x -> plot([[x, x], [x, f(x)], [f(x), f(x)]]):
# Kts. 2.5 Funktioiden yhdistäminen ja iterointi, ss. 64 - 67 (Kuva 2.1 s. 67)
#
# Newton's method:
# ================
# 1-dim:
# ======
N:=x->evalf(x-f(x)/D(f)(x)):
# Käyttöesimerkki:
# f:=x->x*exp(-x)-abs(sin(10*x));
#
# x1:=(N@@5)(0.25);
#
# Operaattorityyli:
 MakeNewtonstep:=proc(flauseke::algebraic,x::name)
    local askel;
    askel:=simplify(x-flauseke/diff(flauseke,x));
    unapply(askel,x);
 end:

# Esim.
#  f:=x-2*sqrt(x);
#  N:=MakeNewtonstep(f,x);                          
#  Iteroidaan, kuten edellä.
#
#Auxiliary function, useful in implementing multidim Newton:
# Multidim Newton can be found on www-page
evm:=proc(mat,symblist,valulist)
   local Mat,i;
   Mat:=mat;
   if type(mat,list) then Mat:=convert(mat,vector) fi;
   for i to nops(symblist) do
       Mat:=subs(symblist[i]=valulist[i],op(Mat))
                           od;
   map(eval,Mat);
   end:

# Esim.
# J:=linalg[jacobian](flist,xlist);Jx0:=evm(J,xlist,x0);

###########
# Ngen -- Generate Matlab code for 1-dim Newton
# Things to add:
#       - include file name in parameter list
#       - Make "terminal" default (or N.m)
# Look at example worksheets in 
#          /p/edu/mat-1.192/k98/luennot-ja-harjoitukset/
#          especially  harj2sol3.mws
# 
###########
Ngen:=proc(f,x)
  local df;
  df:=diff(f,x);
  interface(echo=0);
  system(`rm N.m`);
  writeto(`N.m`);       # Comment for testing
  lprint(`function x1=N(x)`);
  lprint(`f=`,f,`;`);
  lprint(`df=`,df,`;`);
  lprint(`x1 = x -f ./df;`);
  writeto(terminal);
end:

############# 5.4. Numeerinen integrointi #################
# 5.4 s. 117 ->
# Huom: Kirjan kaavassa on virhe jäännöstermissä, tässä korjattu (2. der.)
trapetsi:=proc(f,a,b,n)
  local h,S,R;
  h:=(b-a)/n;
  S:=evalf((h/2)*(f(a)+2*sum(f(a+j*h),j=1..n-1)+f(b)));
  R:=-(h^2*(b-a)/12)*(D@D)(f)(xi);
  [S,R]  # Palutetaan lista: [Intappr, virhe]
end: 
# Esim: 
# f:=x->exp(x)*cos(x);
# trapetsi(f,0,Pi,10);
# 
# Ekonomisempaa saattaa olla tarkastella ensin virhearviota ja vasta
# virhatarkastelun jälkeen laskea integraalin approksimaatio:
trapetsivirhe:=proc(f,a,b,n)
  local h;
  h:=(b-a)/n;
  -(h^2*(b-a)/12)*(D@D)(f)(xi)
end:
#
# Vastaavalla tavalla voidaan toteuttaa Simpsonin kaava.
#
### 5.4. Interpolaatio #########
#
#lag:=(xd,x)->seq(product((x-xd[k])/(xd[j]-xd[k]),k=1..j-1)*
# product((x-xd[k])/(xd[j]-xd[k]),k=j+1..nops(xd)),j=1..nops(xd));
# "Tuotantoversio kirjoitetaan proc-syntaksilla, jotta voidaan lokalisoida
# indeksimuuttujat.
lag:=proc(xd,x)
  local k,j;
  seq(product((x-xd[k])/(xd[j]-xd[k]),k=1..j-1)*
  product((x-xd[k])/(xd[j]-xd[k]),k=j+1..nops(xd)),j=1..nops(xd))
end:
#
# Esim 1: Pisteisiin [a,b,c] liittyvät Lagrangen kertojapolynomit.
# Huom: Ne riippuvat vain "x-datasta".
# lagjono:=lag([a,b,c],x);
# Niiden avulla voidaan kirjoittaa interpolaatiopolynomi, joka datapisteissä
# saa halutut arvot, tässä y[0],y[1], y[2] .
#
# p:=sum(y[i-1]*lagjono[i],i=1..3);
#
# Esim 2:
# lag(linspace(-1,1),x);  # Väli [-1,1] jaetaan 10 osaan ja muod. 
#                         # jakopisteisiin liittyvät lag-polynomit.
#
# HUOM! Maplessa on myös valmis interp, jota kutsutaan:
# interp(xdata,ydata,x);
# 
#  Yllä oleva lag on ennen kaikkea opettelutarkoitukseen (ja lagrangen
#  polynomeja voidaan tarvita muissakin yhteyksissä.
#  Toinen seikka on, että interp-funktioon on rakennettu sisään
#   expand, 
#  mikä ei läheskään aina ole järkevää.
#  lagjono jättää polynomin "Lagrangen muotoon"
#
############# Kompleksianalyysi ###################
# Kompleksitason käyräintegraali
#
 cintf:=proc(f,p,a,b)
   # f on integroitava kompleksifunktio, määriteltävä funktioksi
   # p on polun määrittelevä kompleksifunktio, kuten yllä
   # a ja b ovat parametrivälin alku- ja loppupisteet.
   #
    local intfun,t,dp;
    dp:=diff(p(t),t);
    intfun:=evalc(f(p(t))*dp);
    evalf(Int(intfun,t=a..b));
    end:

############ Tavalliset differntiaaliyhtälöt #####################
#
####################################################################
# Eulerin menetelmä (skalaaritapaus) y'=f(t,y).   y_{n+1}=y_n+h*f(t_n,y_n)
####################################################################
Euler:=proc(f,a,b,ya,m) 
local n,h,t,y;
h:=evalf((b-a)/m);
t[0]:=a;y[0]:=ya;
for n from 0 to m do
  y[n+1]:=y[n]+h*f(t[n],y[n]);
  t[n+1]:=t[n]+h;
od;
[seq([t[n],y[n]],n=0..m)];
end:

#Esim.
###########
#`Example`;
#`f:=(t,y)->t-y^2;`;
#`e3:=Euler(f,0,5,1,3);`;
#`plot(e3);`;

#`For more, see nslib98.mpl-file`;

#Useita piirroksia:
##################
#e6:=Euler(f,0,5,1,6):e12:=Euler(f,0,5,1,12):
#ple3:=plot(e3):ple6:=plot(e6):ple12:=plot(e12):

#with(plots):
#display({ple3,ple6,ple12});
#
############ Vektoriversio ################
#
EulerV:=proc(F,a,b,Ya,m)
local n,h,t,Y;
h := evalf((b - a)/m);
t[0] := a;
Y[0] := Ya;
for n from 0 to m do
        Y[n + 1] := evalm(Y[n] + h*F(t[n],Y[n]));
        t[n + 1] := t[n] + h
    od;
[seq([t[n], op(convert(Y[n],list))], n = 0 .. m)]
end:
#
#
##############################################################
# Parannettu Euler eli Heun (Improved Euler)
##############################################################
Heun:=proc(f,a,b,ya,m) 
local n,h,t,y,ytahti;
h:=evalf((b-a)/m);
t[0]:=a;y[0]:=ya;
for n from 0 to m do
  t[n+1]:=t[n]+h;
  ytahti:=y[n]+h*f(t[n],y[n]);
  y[n+1]:=y[n]+(h/2)*(f(t[n],y[n])+f(t[n+1],ytahti));

od;
[seq([t[n],y[n]],n=0..m)];
end:

############ Vektoriversio ################
#
HeunV:=proc(F,a,b,Ya,m)
 local n, h, t, Y, Ytahti; 
 h := evalf((b-a)/m); 
 t[0] := a; 
 Y[0] := Ya; 
 for n from 0 to m do 
   t[n+1] := t[n]+h; 
   Ytahti := evalm(Y[n]+h*F(t[n],Y[n])); 
   Y[n+1] := evalm(Y[n]+1/2*h*(F(t[n],Y[n])+F(t[n+1],Ytahti)))
 od; 
  [seq([t[n],op(convert(Y[n],list))],n = 0 .. m)] 
 end:

###########################################################
# Runge-Kutta 4
###########################################################
RK4:=proc(f,a,b,ya,m)
local n,h,t,y,k1,k2,k3,k4;
h:=evalf((b-a)/m);
t[0]:=a;y[0]:=ya;
for n from 0 to m do
  k1:=h*f(t[n],y[n]);
  k2:=h*f(t[n]+h/2,y[n]+k1/2);
  k3:=h*f(t[n]+h/2,y[n]+k2/2);
  k4:=h*f(t[n]+h,y[n]+k3);
  y[n+1]:=y[n]+(1/6)*(k1+2*k2+2*k3+k4);
  t[n+1]:=t[n]+h;
od:
[seq([t[n],y[n]],n=0..m)];
end:

##########################################################
# Backward Euler (Implicit Euler)
##########################################################
# BE with Newton iteration, works fine !
##############################################
 BEN:=proc(f,a,b,Y0,N,niter)
 local y0,t0,y1,polyg,G,F,h,i,t1;
 h:=evalf((b-a)/N);
 y0:=Y0;t0:=a;
 polyg:=[a,y0];
 for i from 0 to N do
   t1:=t0+h;
   F:=unapply(y-y0-h*eval(f(t1,y)),y);
   G:=unapply(simplify(y-F(y)/D(F)(y)),y);
   y1:=y0+h*eval(f(t0,y0));    # Initial value for 
                               # iteration (foreuler)
   y1:=(G@@niter)(y1);         # Newton iteration
   polyg:=polyg,[t1,y1];      
   t0:=t1;y0:=y1;
 od:  
 [polyg];
 end:
 
###############################################
# Example:
# f:=(t,y)->-10*(t-1)*y;
# sol:=t->exp(-5*(t-1)^2);
# map(sol,[seq(1+i*.25,i=0..4)]);
# [[BEN(f,1,2,1,4,5)],[Euler(f,1,2,1,4)]];map(matrix,");
#  with(plots):gfor:=plot([Euler(f,1,4,1,4)],color=GREEN):
# gb:=plot([BEN(f,1,4,1,4,5)],color=BLACK):
# gt:=plot(sol,1..5):
# display({gb,gt});
# display(gfor);
####################################################

### Matrix for solving Laplace or Poisson equ using 5-point 
#   difference mesh

lapmatr:=proc(p)
  local B,C,v;
  B:=band([1,-4,1],p);
  C:=diag(seq(B,i=1..p));
  v:=vector(2*p+1,0);
  v[1]:=1;v[2*p+1]:=1;
  evalm(C+band(v,p^2))
end:


## Grid lines for a 2d plot
## Andersson p. 140
##
grid2D:=proc(rx,ry)
  local x,y;
  CURVES(seq([[x,op(1,ry)],[x,op(nops(ry),ry)]],x=rx),
         seq([[op(1,rx),y],[op(nops(rx),rx),y]],y=ry),
         COLOUR(RGB,.5,.5,.5)
        )
end:

## Uniform grid lines, example:
##     grl:=ugrid(0..5,10,-3..7,15)
## 
ugrid:=proc(rx,nx,ry,ny)
  local a,b,c,d,hx,hy,i,j;
  a:=op(1,rx);b:=op(2,rx);
  c:=op(1,ry);d:=op(2,ry);
  hx:=(b-a)/nx;hy:=(d-c)/ny;
  grid2D([seq(a+i*hx,i=0..nx)],[seq(c+j*hy,j=0..ny)])
end:

######################
#Fourier-sarjat
##########################
#
# Jaksollinen laajennus: MapleTech Vol 3 No. 3 1996
#
PeriodicExtender:=proc(f,d::range)
  subs({'F'=f,'L'=lhs(d),'D'=rhs(d)-lhs(d)},
  proc(x::algebraic) local y;
    y:=floor((x-L)/D);
    F(x-y*D);
  end)
end:

# Esim:
#sw:=PeriodicExtender(signum,-1..1);
#plot(sw,-4..4);
#
 parillinenjatko:=f->unapply(piecewise(x>0,f(x),x<0,f(-x)),x):
 paritonjatko:=f->unapply(piecewise(x>0,f(x),x<0,-f(-x)),x):

# Esim. Ei vielä testattu, mutta vrt. työarkkiin fourier.mws
#  fe:=parillinenjatko(x->x);fo:=paritonjatko(x->x);
#  plot([PeriodicExtender(fe,-1..1),PeriodicExtender(fo,-1..1)],-5..5);
#
 parillinen_osa:=f->unapply((1/2)*(f(x)+f(-x)),x):
 pariton_osa:=f->unapply((1/2)*(f(x)-f(-x)),x):

# Parillinen osa on se, josta syntyy cos-sarja, pariton se, josta syntyy
# sin-sarja.
#
#### 7.2. Vektorilaskentaa #########
# Helppokäyttöiset funktiot nuolen piirtämiseen O:sta kohdepisteeseen.
#
snuoli:=u->arrow([0,0],u,0.01,0.05,0.02,color=blue):
pnuoli:=u->arrow([0,0],u,0.01,0.05,0.02,color=red):
#
# Esim: with(plottools):with(plots):
# n1:=snuoli([0,1]):n2:=snuoli([1,0]):n3:=pnuoli([1,1]):
# display([n1,n2,n3]);
#
# Yleiemmin:
nuoli:=(alku,loppu,vari)->arrow(alku,loppu,0.01,0.05,0.02,color=vari):
#
# Esim: display([nuoli([0,0],[0,1],red),nuoli([0,1],[1,1],black),
#               nuoli([1,1],[0,0],blue)]);
#
############ Luku 10 Laplace-muunnokset ###########
#
# Konvoluutio
#
konvoluutio:=(f,g)->unapply(int(f(tau)*g(t-tau),tau=0..t),t):
#
# Esim: Vakiofunktion f(t)=1 ja funktion g(t)=t konvoluutio
#
# h:=konv(t->1,t->t);
#