> 
# Teht. 4
# 
# Muuttujien erottelu
# 
> u:='u':n:='n';  # Varmuuden vuoksi
> osdy:=a^2*diff(u(x,t),x$2)=diff(u(x,t),t):eval(subs(u(x,t)=X(x)*T(t),
> osdy)): #Jaetaan puolittain (a^2*X(x)*T(t)):llä:
> separoitu:=simplify("/(a^2*X(x)*T(t))):xyht:=lhs(separoitu)=-p^2:   #
> Vasen puoli = vakio (negat)
> tyht:=rhs(separoitu)=-p^2:   # Oikea puoli = sama vakio
> xyht:=xyht*X(x):  tyht:=a^2*tyht*T(t): # Kerrotaan nimittäjät pois
> xx:=rhs(dsolve(xyht,X(x))):dsolve(tyht,T(t)): tt:=rhs("):             
>                       
> tt:=subs(_C1=1,tt):   # (Vakio kertaa vakio) yhdistetään yhdeksi
> vakioksi.
>                        # Katso, onko _C1 vai _C2, dsolve voi vaihdella
>                        # järjestystä.
>  # RE1 => ... RE2 => ..., sijoita nämä suoraan (Tässä ei ole mieltä
> "solvata")
> # sin:n 0-kohdat ovat Pi:n monikertoja, kuten tiedetään.
> xx:=subs({_C1=0,_C2=Bn,p=n*Pi/l},xx):
> tt:=subs(p=n*Pi/l,tt);uu:=xx*tt;
# 

                                n := n


                                    2  2   2
                                   a  n  Pi  t
                       tt := exp(- -----------)
                                        2
                                       l


                                           2  2   2
                            n Pi x        a  n  Pi  t
               uu := Bn sin(------) exp(- -----------)
                              l                2
                                              l

# Sinisarjaksi:
# 
> #Bn:=simplify((2/l)*int(F*sin(n*Pi*x/l),x=0..l)); # Tässä F on
> alkuarvofunktio
>                                              # Sijoita se F:n
> paikalle.
> Bn:=simplify((2/l)*(int(50*x*sin(n*Pi*x/l),x=0..1)+int((100-50*x)*sin(
> n*Pi*x/l),x=1..2)));

                                n Pi  /         n Pi \
                          l sin(----) |-1 + cos(----)|
                                 l    \          l   /
               Bn := -200 ----------------------------
                                      2   2
                                     n  Pi

> Bn:=simplify(subs(l=2,Bn));

                       sin(1/2 n Pi) (-1 + cos(1/2 n Pi))
            Bn := -400 ----------------------------------
                                      2   2
                                     n  Pi

> summa:=add(uu,n=1..10);

                              2   2
                 Pi x        a  Pi  t
             sin(----) exp(- --------)
                  l              2
                                l
summa := 400 -------------------------
                          2
                        Pi

                                 2   2
                   Pi x         a  Pi  t
             sin(3 ----) exp(-9 --------)
                    l               2
                                   l
     - 400/9 ----------------------------
                           2
                         Pi

                               2   2
                Pi x          a  Pi  t
          sin(5 ----) exp(-25 --------)
                 l                2
                                 l
     + 16 -----------------------------
                         2
                       Pi

                                2   2
                 Pi x          a  Pi  t
           sin(7 ----) exp(-49 --------)
                  l                2
       400                        l
     - --- -----------------------------
       49                 2
                        Pi

                                2   2
                 Pi x          a  Pi  t
           sin(9 ----) exp(-81 --------)
                  l                2
       400                        l
     + --- -----------------------------
       81                 2
                        Pi

> summa:=subs(a=sqrt(2.9*10^(-5)),l=2,summa);

                                              -5   2
             sin(1/2 Pi x) exp(-.7250000000 10   Pi  t)
summa := 400 ------------------------------------------
                                  2
                                Pi

                                                  2
             sin(3/2 Pi x) exp(-.00006525000000 Pi  t)
     - 400/9 -----------------------------------------
                                  2
                                Pi

                                              2
          sin(5/2 Pi x) exp(-.0001812500000 Pi  t)
     + 16 ----------------------------------------
                              2
                            Pi

                                               2
       400 sin(7/2 Pi x) exp(-.0003552500000 Pi  t)
     - --- ----------------------------------------
       49                      2
                             Pi

                                               2
       400 sin(9/2 Pi x) exp(-.0005872500000 Pi  t)
     + --- ----------------------------------------
       81                      2
                             Pi

> 
> plot3d(summa,x=0..2,t=0..5000,axes=boxed);

> ufun:=unapply(summa,x,t);

                                                       -5   2
                      sin(1/2 Pi x) exp(-.7250000000 10   Pi  t)
ufun := (x, t) -> 400 ------------------------------------------
                                           2
                                         Pi

                                                  2
             sin(3/2 Pi x) exp(-.00006525000000 Pi  t)
     - 400/9 -----------------------------------------
                                  2
                                Pi

                                              2
          sin(5/2 Pi x) exp(-.0001812500000 Pi  t)
     + 16 ----------------------------------------
                              2
                            Pi

                                               2
       400 sin(7/2 Pi x) exp(-.0003552500000 Pi  t)
     - --- ----------------------------------------
       49                      2
                             Pi

                                               2
       400 sin(9/2 Pi x) exp(-.0005872500000 Pi  t)
     + --- ----------------------------------------
       81                      2
                             Pi

> plot({seq(ufun(x,t),t=[0,1000,3600,10000])},x=0..2);
> unapply(x^2,x);


> 
> 

>