Lämpöyhtälö muuttujat erotellen

Numsym 15.2.2005

Alustukset

>    restart:

>    with(plots): setoptions3d(orientation=[-120,50],axes=box):

Warning, the name changecoords has been redefined

Esim

KRE s. 603 Exa 1

>    f:=x->100*sin(Pi*x/80); plot(f,0..80);

f := proc (x) options operator, arrow; 100*sin(1/80*Pi*x) end proc

[Maple Plot]

>    u:=(x,t)->100*sin(Pi*x/80)*exp(-0.001785*t);

u := proc (x, t) options operator, arrow; 100*sin(1/80*Pi*x)*exp(-.1785e-2*t) end proc

>   

>    plot3d(u(x,t),x=0..80,t=0..400,axes=box);

[Maple Plot]

 "Aikasiivuparvi" :

>    [seq(u(x,t),t=[seq(i*0.1,i=0..5)])];

[100.*sin(1/80*Pi*x), 99.98215159*sin(1/80*Pi*x), 99.96430637*sin(1/80*Pi*x), 99.94646434*sin(1/80*Pi*x), 99.92862548*sin(1/80*Pi*x), 99.91078982*sin(1/80*Pi*x)]

Tällainen lausekkeiden lista voidaan antaa suoraan plot: lle, kuten tiedämme: Otetaan samantien enemmän aikaviipaleita.

>    plot([seq(u(x,t),t=[seq(i*30,i=0..10)])],x=0..80);

[Maple Plot]

>    animate(u(x,t),x=0..80,t=0..400,frames=30);

[Maple Plot]

Exa 2. Muuten sama, mutta nyt

>    restart:with(plots): setoptions3d(orientation=[-120,50],axes=box):

Warning, the name changecoords has been redefined

>    f:=x->100*sin(3*Pi*x/80); plot(f,0..80);

f := proc (x) options operator, arrow; 100*sin(3/80*Pi*x) end proc

[Maple Plot]

>    lambda[3]:=c*3*Pi/L;

lambda[3] := 3*c*Pi/L

>    c:=sqrt(1.158):L:=80: lambda[3];

.4035390315e-1*Pi

>    u:=(x,t)->100*sin(3*Pi*x/L)*exp(-lambda[3]^2*t);

u := proc (x, t) options operator, arrow; 100*sin(3*Pi*x/L)*exp(-lambda[3]^2*t) end proc

>    plot3d(u(x,t),x=0..80,t=0..200,axes=box);

[Maple Plot]

 "Aikasiivuparvi" :

>    #[seq(u(x,t),t=[seq(i*0.1,i=0..5)])];

Tällainen lausekkeiden lista voidaan antaa suoraan plot: lle, kuten tiedämme: Otetaan samantien enemmän aikaviipaleita.

>    plot([seq(u(x,t),t=[seq(i*20,i=0..10)])],x=0..80);

[Maple Plot]

>    animate(u(x,t),x=0..80,t=0..200,frames=30);

[Maple Plot]

>   

>   

Muuttujien erottelu

>    restart:  
with(LinearAlgebra): with(plots):

Warning, the name changecoords has been redefined

>    setoptions3d(orientation=[-120,50],axes=box):

>    #read("c:\\usr\\heikki\\numsym05\\maple\\ns05.mpl");

>    read("e:\\ns05.mpl");

>    # read("/p/edu/mat-1.192/ns05.mpl");

>    osdy:=c^2*diff(u(x,t),x$2)=diff(u(x,t),t);eval(subs(u(x,t)=F(x)*G(t), osdy));

osdy := c^2*diff(u(x,t),`$`(x,2)) = diff(u(x,t),t)

c^2*diff(F(x),`$`(x,2))*G(t) = F(x)*diff(G(t),t)

>    separoitu:=simplify(%/(c^2*F(x)*G(t)));

separoitu := 1/F(x)*diff(F(x),`$`(x,2)) = 1/c^2/G(t)*diff(G(t),t)

>   

>    xyht:=lhs(separoitu)=-p^2;
tyht:=rhs(separoitu)=-p^2;

xyht := 1/F(x)*diff(F(x),`$`(x,2)) = -p^2

tyht := 1/c^2/G(t)*diff(G(t),t) = -p^2

>    FDY:=lhs(xyht)*F(x)=rhs(xyht)*F(x);

FDY := diff(F(x),`$`(x,2)) = -p^2*F(x)

>    GDY:=c^2*lhs(tyht)*G(t)=c^2*rhs(tyht)*G(t);

GDY := diff(G(t),t) = -c^2*p^2*G(t)

>   

>   

>    ff:=rhs(dsolve(FDY,F(x)));
dsolve(GDY,G(t)): gg:=rhs(%);

ff := _C1*sin(p*x)+_C2*cos(p*x)

gg := _C1*exp(-c^2*p^2*t)

>    ff:=subs(_C2=A,_C1=B,ff); gg:=subs(_C1=1,gg);

ff := B*sin(p*x)+A*cos(p*x)

gg := exp(-c^2*p^2*t)

RE1 & RE2 ==>

>    ff:=subs(A=0,B=1,p=n*Pi/L,ff);

ff := sin(n*Pi/L*x)

>    gg:=subs(p=n*Pi/L,gg);

gg := exp(-c^2*n^2*Pi^2/L^2*t)

>    uu:=ff*gg;

uu := sin(n*Pi/L*x)*exp(-c^2*n^2*Pi^2/L^2*t)

>    u:=sum(B[n]*ff*gg,n=1..infinity);

u := sum(B[n]*sin(n*Pi/L*x)*exp(-c^2*n^2*Pi^2/L^2*t),n = 1 .. infinity)

Sinisarjaksi:

Alkuehto saadaan toteutumaan valitsemalla B[n]:t AE-funktion f  (2L-jaksoisen) sinisarjan kertoimiksi.

Määritellään sinisarjan kertoimet ja annetaan alkuarvofunktion f sekä jakson puolikkaan L olla

vapaita, myöhemmin annettavia symboleja (vrt. fourier.mws-työarkki).

>    f:='f':L:='L':n:='n':Bn:=(2/L)*Int(f(x)*sin(n*Pi*x/L),x=0..L);              

Bn := 2/L*Int(f(x)*sin(n*Pi/L*x),x = 0 .. L)

Otetaan esimerkiksi KRE s. 648, Exa 3

>    f:=x->piecewise(x<L/2,x,x>L/2,L-x);L:=Pi;

f := proc (x) options operator, arrow; piecewise(x < 1/2*L,x,1/2*L < x,L-x) end proc

L := Pi

>    plot(f,0..Pi);

[Maple Plot]

pi :n pituinen sauva sopiii matemaatikoille.

>   

>    bn:=simplify(value(Bn));

bn := -2/Pi*(sin(Pi*n)-2*sin(1/2*Pi*n))/n^2

>    bn:=trigsiev(%,n);

bn := 4/Pi*sin(1/2*Pi*n)/n^2

>    uu;

sin(n*x)*exp(-c^2*n^2*t)

>    un:=bn*uu;

un := 4/Pi*sin(1/2*Pi*n)/n^2*sin(n*x)*exp(-c^2*n^2*t)

>    summaf:=(x,t,N)->add(un,n=1..N);

summaf := proc (x, t, N) options operator, arrow; add(un,n = 1 .. N) end proc

>    summaf(x,t,4);

4/Pi*sin(x)*exp(-c^2*t)-4/9*1/Pi*sin(3*x)*exp(-9*c^2*t)

>    summaf(x,t,10);

4/Pi*sin(x)*exp(-c^2*t)-4/9*1/Pi*sin(3*x)*exp(-9*c^2*t)+4/25*1/Pi*sin(5*x)*exp(-25*c^2*t)-4/49*1/Pi*sin(7*x)*exp(-49*c^2*t)+4/81*1/Pi*sin(9*x)*exp(-81*c^2*t)

>    c:=1:

>    plot3d(summaf(x,t,5),x=0..Pi,t=0..2,axes=box);

[Maple Plot]

Summafunktion "aikasiivuparvi" :

>    [seq(summaf(x,t,5),t=[seq(i*0.1,i=0..5)])];

[4./Pi*sin(x)-.4444444444/Pi*sin(3*x)+.1600000000/Pi*sin(5*x), 3.619349672/Pi*sin(x)-.1806976265/Pi*sin(3*x)+.1313359978e-1/Pi*sin(5*x), 3.274923012/Pi*sin(x)-.7346617252e-1/Pi*sin(3*x)+.1078071520e-2/...
[4./Pi*sin(x)-.4444444444/Pi*sin(3*x)+.1600000000/Pi*sin(5*x), 3.619349672/Pi*sin(x)-.1806976265/Pi*sin(3*x)+.1313359978e-1/Pi*sin(5*x), 3.274923012/Pi*sin(x)-.7346617252e-1/Pi*sin(3*x)+.1078071520e-2/...
[4./Pi*sin(x)-.4444444444/Pi*sin(3*x)+.1600000000/Pi*sin(5*x), 3.619349672/Pi*sin(x)-.1806976265/Pi*sin(3*x)+.1313359978e-1/Pi*sin(5*x), 3.274923012/Pi*sin(x)-.7346617252e-1/Pi*sin(3*x)+.1078071520e-2/...
[4./Pi*sin(x)-.4444444444/Pi*sin(3*x)+.1600000000/Pi*sin(5*x), 3.619349672/Pi*sin(x)-.1806976265/Pi*sin(3*x)+.1313359978e-1/Pi*sin(5*x), 3.274923012/Pi*sin(x)-.7346617252e-1/Pi*sin(3*x)+.1078071520e-2/...
[4./Pi*sin(x)-.4444444444/Pi*sin(3*x)+.1600000000/Pi*sin(5*x), 3.619349672/Pi*sin(x)-.1806976265/Pi*sin(3*x)+.1313359978e-1/Pi*sin(5*x), 3.274923012/Pi*sin(x)-.7346617252e-1/Pi*sin(3*x)+.1078071520e-2/...

Tällainen lausekkeiden lista voidaan antaa suoraan plot: lle, kuten tiedämme: Otetaan samantien enemmän aikaviipaleita.

>    plot([seq(summaf(x,t,40),t=[seq(i*0.1,i=0..10)])],x=0..Pi);

[Maple Plot]

>    animate(summaf(x,t,40),x=0..Pi,t=0..10,frames=30);

[Maple Plot]

>   

Eristetyt päät, adiabaattiset reunaehdot

>    restart:  
with(plots):setoptions3d(orientation=[-120,50],axes=box):

Warning, the name changecoords has been redefined

>    #read("c:\\usr\\heikki\\numsym05\\maple\\ns05.mpl");

>    read("e:\\ns05.mpl");

>    uu:=(x,t,n)->cos(n*Pi*x/L)*exp(-(c*n*Pi/L)^2*t);

uu := proc (x, t, n) options operator, arrow; cos(n*Pi*x/L)*exp(-c^2*n^2*Pi^2/L^2*t) end proc

>    usum:=(x,t,N)->add(a(n)*uu(x,t,n),n=0..N);

usum := proc (x, t, N) options operator, arrow; add(a(n)*uu(x,t,n),n = 0 .. N) end proc

>    usum(x,t,3);

a(0)+a(1)*cos(Pi*x/L)*exp(-c^2*Pi^2/L^2*t)+a(2)*cos(2*Pi*x/L)*exp(-4*c^2*Pi^2/L^2*t)+a(3)*cos(3*Pi*x/L)*exp(-9*c^2*Pi^2/L^2*t)

>    f:='f':L:='L':n:='n':an:=(2/L)*Int(f(x)*cos(n*Pi*x/L),x=0..L);              

an := 2/L*Int(f(x)*cos(n*Pi*x/L),x = 0 .. L)

Esimerkiksi KRE s. 64x, Exa 4

>    f:=x->piecewise(x<L/2,x,x>L/2,L-x);L:=Pi;

f := proc (x) options operator, arrow; piecewise(x < 1/2*L,x,1/2*L < x,L-x) end proc

L := Pi

>    plot(f,0..Pi);

[Maple Plot]

>   

>    a:=n->2/L*int(f(x)*cos(n*Pi*x/L),x=0..L);

a := proc (n) options operator, arrow; 2/L*int(f(x)*cos(n*Pi*x/L),x = 0 .. L) end proc

>    a(0):=1/L*int(f(x),x=0..L);

a(0) := 1/4*Pi

>    a(n);

2/Pi*(-cos(Pi*n)+2*cos(1/2*Pi*n)-1)/n^2

>    a(0);

1/4*Pi

>    seq(a(n),n=0..10);

1/4*Pi, 0, -2/Pi, 0, 0, 0, -2/9/Pi, 0, 0, 0, -2/25/Pi

>    a(1);

0

>    usum(x,t,5);

1/4*Pi-2/Pi*cos(2*x)*exp(-4*c^2*t)

>    c:=1:

>    plot3d(usum(x,t,20),x=0..Pi,t=0..2,axes=box);

[Maple Plot]

Summafunktion "aikasiivuparvi" :

>    #[seq(usum(x,t,20),t=[seq(i*0.1,i=0..5)])];

Tällainen lausekkeiden lista voidaan antaa suoraan plot: lle, kuten tiedämme: Otetaan samantien enemmän aikaviipaleita.

>    plot([seq(usum(x,t,10),t=[seq(i*0.1,i=0..10)])],x=0..Pi);

[Maple Plot]

>    animate(usum(x,t,20),x=0..Pi,t=0..2,frames=30);

[Maple Plot]

>   

>