### maple/v301.mpl    25.9.2001
# Tässä on joitakin viime kevään työkaluja, mutta ei kaikkia,
# loput:  maple/v201.mpl
#
# read("c:\\windows\\mymaple\\v301.mpl");  (PC-WIN:ssä)  
# read("/p/edu/mat-1.414/maple/v301.mpl");
# read("/home/apiola/opetus/peruskurssi/v2-3/201/maple/v301.mpl");  HA:n priv.

############################################
##Datan muokkausta, yleisiä apufunktioita: 
############################################

v2l:=vek->convert(vek,list):    # HAM:ssa painovirhe (vec)
l2v:=lis->convert(lis,vector):

#Liitetään kaksi vektoria peräkkäin: 

vappend:=(v1,v2)->vector([op(convert(v1,list)),op(convert(v2,list))]):

#Matriisin jonouttaminen pitkäksi vektoriksi 

m2l:=M->map(op,convert(M,listlist)):  # Vastaa Matlabin M(:), mutta tässä
                                      # riveittäin 

#Diskretointi vrt. Matlab-tuttavamme 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 [HAM] s. 117. Vaihdoimme oletusarvon 10:een,
#Matlab:n 100:n sijasta. 

         ###################
         #Konformikuvauksia#
         ###################
gympyra:=(z0,r)->complexplot(z0+r*exp(I*Theta),Theta=0..2*Pi):
gympyrakuva:=(z0,r,f)->complexplot(f(z0+r*exp(I*Theta)),Theta=0..2*Pi):

#display([seq(gympyra(I+0.1*k,1),k=0..20)],insequence=true);
#display([seq(gympyra(0,0.1*k),k=0..10)]);

#gsade:=(Theta,r1,r2)->complexplot(r*exp(I*Theta),r=r1..r2):
#gsadekuva:=(Theta,r1,r2,f)->complexplot(f(r*exp(I*Theta)),r=r1..r2):
#display([seq(gsade(k*2*Pi/10,0.5,3),k=0..9)]);

#nymp:=7: nsateet:=15: rmax:=4: hr:=rmax/nymp: hTheta:=2*Pi/nsateet:
#display([seq(gsade(k*hTheta,0,rmax),k=0..nsateet),seq(gympyra(0,k*hr),
#k=0..nymp)]);

       ##############
       # Iteraatiot #
       ##############

porras:=x->plot([[x,x],[x,g(x)],[g(x),g(x)]]):
newtonkulma:=x->plot([[x,0],[x,f(x)],[N(x),0]],color=black):
NNN:=x->x-f(x)/D(f)(x):   # Pelkkä N on liian salakavala!
# N:=NNN                   


       #################
       # Interpolaatio #
       #################

# Lagrangen kertojapolynomi:
L:=proc(j,xd,x)
local oso,nimi,i,j1;
j1:=j+1;
oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j1]);
nimi:=subs(x=xd[j1],oso);
oso/nimi;
end:

#Parametrit 
#  j  -- antaa L:n indeksin, eli kyseessä on Lj
#        j=0..n, missä n on xd-listan pituus - 1 
#  xd -- xdata, [x0,x1,...,xn], lista (n+1 pistettä)
#  x  -- polynomin muuttuja

#Lagrangen polynomijono muodostetaan nyt näin, kun data on listassa xd

#    seq(L(j,xd,x),j=0..nops(xd)-1);

#################################################
# Interpolaatiopolynomi Lagrangen menetelmällä:
#################################################
lagint:=(xd,yd,x)->add(yd[j+1]*L(j,xd,x),j=0..nops(xd)-1):   

# lagint([a,b,c],[ya,yb,yc],x);
# xd:=[a,b,c]: lagint(xd,map(f,xd),x);

##########################################################################
# interplot -- piirretään datapisteisiin liittyvä interpolaatiopolynomi
# ja datapisteet samaan kuvaan
# Tämä on niin usein esiintyvä tehtävä, että kannatti kirjoittaa proc:ksi.
#  HA 7.2.01
#######################################################################

 interplot:=proc(xd,yd)
 local xyd,p,x;
 xyd:=zip((x,y)->[x,y],xd,yd);
 p:=interp(xd,yd,x);
 plots[display]([plot(p,x=xd[1]..xd[nops(xd)],color=blue),plot(xyd,style=point,symbol=cross,symbolsize=20,color=black)]);
 end:

# Esim: xd:=[...]: f:=x->...: yd:=map(f,xd):
#       interplot(xd,yd);
#       display([interplot(xd,yd),f(x),x=xd[1]..xd[nops(xd)]]);

######################################################################
#  Pelkkä datan piirto, tätä tarvitaan usein. (Tässä on riisuttu interplotista
#  muut osat pois )
#  Data annetaan kahtena listana; xd ja yd
#
 dataplot:=proc(xd,yd)
 local xyd;
 xyd:=zip((x,y)->[x,y],xd,yd);
 plot(xyd,style=point,symbol=cross,symbolsize=20,color=black);
 end:

######################################################
#
# 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]);
#
# Yleisemmin:
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)]);
#


       #########################
       # Grafiikkafunktioita #
       #########################

## Koordinaattiristikko 19.4.01 

vaakajana:=(x1,x2,y)->plot([[x1,y],[x2,y]],thickness=1,color=grey):
pystyjana:=(y1,y2,x)->plot([[x,y1],[x,y2]],thickness=1,color=grey):

hila:=(a,b,m,c,d,n)->display([seq(vaakajana(a,b,c+_k*(d-c)/n),_k=0..n),seq(pystyjana(c,d,a+_kk*(b-a)/m),_kk=0..m)]):


#Käytämme seq-muuttujia _k, _kk,koska -> tyylissä ei voi lokalisoida
# (tietääkseni)

# Esim:
#      hila(-1,1,20,3,10,4);       
#      (siistiä!)

#########################################################################
# V3-funktioita syksy 2001 
################################
# HAM-kirjan koodeja: 
#
# Picard-Lindelöfin iteraatio (Alkuarvotehtävän olemassolo- ja 1-käs.)
picklind1:=(fxt,t,x,xfun,t0,x0)->
          x0+int(subs(x=xfun,t=tau,fxt),tau=t0..t):

PL:=xfun->picklind1(x^2+t^2,t,x,xfun,0,1);

#(PL@@3)(1); 
#seq((PL@@k)(1),k=1..4);
#plot([%],t=0..1);


#################################################
# Eulerin menmetelmä diffyhtälöön (AA)
#
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:

#Example
###########
#`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});
###########################################################
# 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);
####################################################
#
# Värähtelevä jousi 
# Lähde: Matta-projekti: http://matta.hut.fi/matta/deltampl/delta1/index.html
#           -> sovellukset -> värähtelevä jousi.
#
# Tässä on yo. sivulta otettu koodi. Alla tulevat määritellyksi funktiot
#     jousi, kappale, varahtelija
######### Istunto, poista kommenteja, jos haluat suorittaa .#############
# restart; 
#yht:= m*diff(x(t), t$2)+k*x(t)=0; 
#alkuehdot:= x(0)=0, D(x)(0)=10; 
#rtk:= dsolve({yht, alkuehdot}, x(t)); 
#rata1:=subs({k=1, m=1}, rtk); 
#plot(rhs(rata1), t=0..4*Pi); 

#with(plots): with(plottools): 

jousi:= proc(ala, yla, paa, jaksot)
  local ampl, jakso;
  ampl:= 3; 
  jakso:= (yla-ala-2*paa)/jaksot;
  plot(
  [[ampl*sin(2*Pi*r/jakso),
  r+ala+paa, r=0..jaksot*jakso],
  [0, ala+r, r=0..paa],
  [0, yla-r, r=0..paa]], color=black);
  end:

  kappale:= proc(hor, lev, kor)
  rectangle([-0.5*lev, hor+0.5*kor],
  [0.5*lev, hor-0.5*kor], color=gray)
  end:

  varahtelija:= proc(y0, t0)
  local lev, kor, paa, lepo, y1;
  lev:= 20; kor:= 10; paa:= 5; lepo:= 50;
  y1:= -lepo+subs(t=t0/(2*Pi), subs(y0,x(t)));
  display(
  rectangle([-50, 0], [50, 5], color=gray),
  kappale(y1, lev, kor), 
  jousi(y1+kor/2, 0, 5, 5));
  end: 

#display(seq(varahtelija(rata1, u), u=0..39), insequence=true, axes=none,
#  scaling=constrained);