A házi feladat dokumentációja:

Adatok:

Numerikus- és szimbolikus számítások házi feladat

Név: Stippinger Marcell

BME-TTK mérnök-fizikus szak, 1. éves hallgató

Neptun kód: XXXXXX

A választott feladat:

A fizikában gyakori eset, hogy a törvényeket "pontosan" ismerjük, fel tudjuk írni a differenciálegyenleteket, de a feladat megoldásához ezzel nem jutunk közelebb.

Az egyenlet által leírt mozgás lehet kaotikus, azaz a kezdeti feltételek kis megváltoztatásaesetén teljesen más eredményt kapunk. Erre példaként a Föld időjárását szokták említeni. De nem kell olyan nagy léptékben gondolkozni, bárki készíthet otthon is egy kaotikus ingát (csillapodás esetén nem a végállás, hanem a befutott út tekintetében alkalmazhatjuk rá, hogy kaotikus). Ez a tanulmány ennek modellezéséről szól.

Ez a maple worksheet egy kaotikus inga mozgását modellezi le. Az inga lényege:

Egy asztalon erős mágneseket rögzítünk, melyek hatnak a fölöttük mozgó mágnesből készült matematikainak tekinthető ingára. Paraméterként módosítható a mágnesek száma, helyzete, hossza, erőssége, az inga kezdeti pozíciója és sebessége.

A mágnességtan felhasznált törvényei:

A rúdmágnes felfogható mágneses dipólusként, melynek hossza a mágnes hosszának 5/6-a, a pólusok pedig a mágnes hosszának 1/12 résszével beljebb vannak a mágnes végétől:

d[m] = p*5*l/6

Az adott ponban érvényes mágneses indukcióvektor a mágneses Coulomb törvény és a szuperpozícó elve alapján:

B = [1/4*Pi]*sum(p[i]*r/abs(r)/r^2, i = (1 .. n))

Adott H térben a p erősségű pólusra ható erő:

(F = p*H) = p*B/mu

Az inga mozgásegyenlete (p pozitív érték):

m*diff(r[x](t), t, t)+k*diff(r[x](t), t)+m*g/l*r[x](t) = {p*B(r[pozitiv_polus])/mu-p*B(r[negativ_polus])/mu}[x]

és

m*diff(r[y](t), t, t)+k*diff(r[y](t), t)+m*g/l*r[y](t) = {p*B(r[pozitiv_polus])/mu-p*B(r[negativ_polus])/mu}[y]

Ahhoz, hogy a feladat megoldása kicsit egyszerűbb, és az ábrák szebbek legyenek, feltesszük, hogy az inga a lengés során nem távolodik el az asztaltól (vagy az aztal hozzá van görbítve az ingához).

Felhasznált irodalom:

Budó Ágoston: Kísérleti fizika II.

Középiskolai fizikatanulmányok: iskolai rúdmágnes póluserőssége.

A maple inicializálása a feladathoz

> restart;

A feladaban felhasznált számítási eszközök:

A vektorok kezeléséhez a VectorCalculus Package-t választottam, mert kb. 5-ször gyorsabb, mint a LinearAlgebra Package. A differenciálegyenlet numerikus megoldását és ábrázolását a DEtools Package eszközeiel végzem.

> with(VectorCalculus): with(DEtools):

Warning, the assigned names `<,>` and `<|>` now have a global binding

Warning, these protected names have been redefined and unprotected: `*`, `+`, `.`, D, Vector, diff, int, limit, series

1. beállítás: a két szélső mágnes vonzó, a középsők taszítóak

A mágnesekre vonatkozó konstansok

A mágnességtan fikikai konstansai:

> mu0:=4*Pi/10^7;   # Vs/Am
mur:=1;           # 1

mu:=mu0*mur;      # Vs/Am

mu0 := 1/2500000*Pi

mur := 1

mu := 1/2500000*Pi

A szimulációban használt mágnesek száma (az ingával együtt):

> magn_db:=5;

magn_db := 5

A mágnesek póluserőssége (pozitív számként adjuk meg ezeket, az első az inga):

> magn_p:=
[6.6E-5,

3.6E-5,

3.6E-5,

3.6E-5,

3.6E-5];

# Weber = Wb = V*s

magn_p := [0.66e-4, 0.36e-4, 0.36e-4, 0.36e-4, 0.36e-4]

A mágnesek hossza (a pozitív hossz jelenti, hogy fölfelé áll a pozitív pólus, az első az inga):

> magn_l:=
[+0.06,

+0.06,

-0.06,

-0.06,

+0.06];

# meter = m

magn_l := [0.6e-1, 0.6e-1, -0.6e-1, -0.6e-1, 0.6e-1]

A mágnesek középpontjainak koordinátái (az első az inga):

> magn_choords:=
[Vector([   0.0 ,   0.0 ,    0.012]),

Vector([   0.15,   0.05,   -0.050]),

Vector([   0.05,  -0.05,   -0.050]),

Vector([  -0.05,   0.05,   -0.050]),

Vector([  -0.15,  -0.05,   -0.050])];

# meter = m

magn_choords := [Vector[column]([[0.], [0.], [0.12e-1]], [

Mágnességtani számítások

Az asztal fölötti adott pontban érvényes térerősséget B eljárás számolja ki:

a for ciklus összegzi az asztalhoz rögzített mágnesek hatását B0 eredővektorba,

r1 az i. mágnes pozitív pólusának, r2 a negatív pólusának helyvektora.

> B:=proc(r::Vector)
local B0,B1,B2,r1,r2,i;

B0:=Vector(3,0);

for i from 2 to magn_db do

 r1:=r-magn_choords[i];

 r1[3]:=r1[3]-5/12*magn_l[i];

 B1:=(  magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r1,r1))^3) )  )*r1;

 r2:=r-magn_choords[i];

 r2[3]:=r2[3]+5/12*magn_l[i];

 B2:=-( magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r2,r2))^3) )  )*r2;

 B0:=B0+B1+B2;

od;

return B0;

end proc:

Grafikus megjelenítés kiegészítései

A grafikonok színezéséhez itt egyedi függvényt lehet definiálni.

> szin:=proc(a)
return a;

end proc;

szin := proc (a) return a end proc

A feladat fizikai megfogalmazása

A mozgás differenciálegyenletei:

diff(x(t), t, t)*m+diff(x(t), t)*k+x(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[1]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))*magn_p[1])/mu

és

diff(y(t), t, t)*m+diff(y(t), t)*k+y(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[2]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))[2]*magn_p[1])/mu

ahol m az inga tömege, l az inga hossza, g a gravivációs gyorsulási együttható, k a közegellenállás.

A megoldások ábrázolása

> m:=0.2; k:=.01; l:=5; g:=9.81; z:=magn_choords[1][3];
me1:=

 diff(x(t),t,t)*m+diff(x(t),t)*k+x(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[1]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[1]*magn_p[1] )/mu:

me2:=

 diff(y(t),t,t)*m+diff(y(t),t)*k+y(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[2]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[2]*magn_p[1] )/mu:

m := .2

k := 0.1e-1

l := 5

g := 9.81

z := 0.12e-1

> x0:=0.2: y0:=0.0: vx:=-0.1: vy:=0.01: tm:=30: xm:=0.3: ym:=0.3:
a1:=DEplot3d([me1,me2],

      [x(t),y(t)],t=0..tm,x=-xm..xm,y=-ym..ym,

      [[x(0)=x0,y(0)=y0,(D)(x)(0)=vx,(D)(y)(0)=vy]],

      stepsize=0.01,linecolor=szin(t),thickness=2):

list_polys:=[seq([[0 ,magn_choords[i][1],magn_choords[i][2] ],

                 [tm,magn_choords[i][1],magn_choords[i][2] ]],

                i=2..magn_db)]:

a4:=plots[polygonplot3d](list_polys, thickness=3):

plots[display]({a1,a4},orientation=[90,180]);

plots[display]({a1,a4},orientation=[-90,90]);

plots[display]({a1,a4},orientation=[1,88]);

[Plot]

[Plot]

[Plot]

>

Érdekesség, hogy mind a vonzó, mind a taszító mágnest kb. ugyanannyira közelíti meg az inga:

a vonzónál (fel)gyorsulva halad el - ezért megy el mellette, a taszító ellöki magától.

2. beállítás: mind a négy mágnes taszító

A mágnesekre vonatkozó konstansok

A mágnességtan fikikai konstansai:

> mu0:=4*Pi/10^7;   # Vs/Am
mur:=1;           # 1

mu:=mu0*mur;      # Vs/Am

mu0 := 1/2500000*Pi

mur := 1

mu := 1/2500000*Pi

A szimulációban használt mágnesek száma (az ingával együtt):

> magn_db:=5;

magn_db := 5

A mágnesek póluserőssége (pozitív számként adjuk meg ezeket, az első az inga):

> magn_p:=
[6.6E-5,

3.6E-5,

3.6E-5,

3.6E-5,

3.6E-5];

# Weber = Wb = V*s

magn_p := [0.66e-4, 0.36e-4, 0.36e-4, 0.36e-4, 0.36e-4]

A mágnesek hossza (a pozitív hossz jelenti, hogy fölfelé áll a pozitív pólus, az első az inga):

> magn_l:=
[+0.06,

-0.06,

-0.06,

-0.06,

-0.06];

# meter = m

magn_l := [0.6e-1, -0.6e-1, -0.6e-1, -0.6e-1, -0.6e-1]

A mágnesek középpontjainak koordinátái (az első az inga):

> magn_choords:=
[Vector([   0.0 ,   0.0 ,    0.012]),

Vector([   0.15,   0.05,   -0.050]),

Vector([   0.05,  -0.05,   -0.050]),

Vector([  -0.05,   0.05,   -0.050]),

Vector([  -0.15,  -0.05,   -0.050])];

# meter = m

magn_choords := [Vector[column]([[0.], [0.], [0.12e-1]], [

Mágnességtani számítások

Az asztal fölötti adott pontban érvényes térerősséget B eljárás számolja ki:

a for ciklus összegzi az asztalhoz rögzített mágnesek hatását B0 eredővektorba,

r1 az i. mágnes pozitív pólusának, r2 a negatív pólusának helyvektora.

> B:=proc(r::Vector)
local B0,B1,B2,r1,r2,i;

B0:=Vector(3,0);

for i from 2 to magn_db do

 r1:=r-magn_choords[i];

 r1[3]:=r1[3]-5/12*magn_l[i];

 B1:=(  magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r1,r1))^3) )  )*r1;

 r2:=r-magn_choords[i];

 r2[3]:=r2[3]+5/12*magn_l[i];

 B2:=-( magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r2,r2))^3) )  )*r2;

 B0:=B0+B1+B2;

od;

return B0;

end proc:

Grafikus megjelenítés kiegészítései

A grafikonok színezéséhez itt egyedi függvényt lehet definiálni.

> szin:=proc(a)
return a;

end proc;

szin := proc (a) return a end proc

A feladat fizikai megfogalmazása

A mozgás differenciálegyenletei:

diff(x(t), t, t)*m+diff(x(t), t)*k+x(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[1]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))*magn_p[1])/mu

és

diff(y(t), t, t)*m+diff(y(t), t)*k+y(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[2]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))[2]*magn_p[1])/mu

ahol m az inga tömege, l az inga hossza, g a gravivációs gyorsulási együttható, k a közegellenállás.

A megoldások ábrázolása

> m:=0.2; k:=.01; l:=5; g:=9.81; z:=magn_choords[1][3];
me1:=

 diff(x(t),t,t)*m+diff(x(t),t)*k+x(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[1]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[1]*magn_p[1] )/mu:

me2:=

 diff(y(t),t,t)*m+diff(y(t),t)*k+y(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[2]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[2]*magn_p[1] )/mu:

m := .2

k := 0.1e-1

l := 5

g := 9.81

z := 0.12e-1

> x0:=0.2: y0:=0.0: vx:=-0.1: vy:=0.01: tm:=30: xm:=0.3: ym:=0.3:
a1:=DEplot3d([me1,me2],

      [x(t),y(t)],t=0..tm,x=-xm..xm,y=-ym..ym,

      [[x(0)=x0,y(0)=y0,(D)(x)(0)=vx,(D)(y)(0)=vy]],

      stepsize=0.01,linecolor=szin(t),thickness=2):

list_polys:=[seq([[0 ,magn_choords[i][1],magn_choords[i][2] ],

                 [tm,magn_choords[i][1],magn_choords[i][2] ]],

                i=2..magn_db)]:

a4:=plots[polygonplot3d](list_polys, thickness=3):

plots[display]({a1,a4},orientation=[90,180]);

plots[display]({a1,a4},orientation=[-90,90]);

plots[display]({a1,a4},orientation=[1,88]);

[Plot]

[Plot]

[Plot]

>

A kezdeti feltételek érdekes, változatos pályát határoznak meg.

3. beállítás: négy, egy egyenesbe eső mágnes

A mágnesekre vonatkozó konstansok

A mágnességtan fikikai konstansai:

> mu0:=4*Pi/10^7;   # Vs/Am
mur:=1;           # 1

mu:=mu0*mur;      # Vs/Am

mu0 := 1/2500000*Pi

mur := 1

mu := 1/2500000*Pi

A szimulációban használt mágnesek száma (az ingával együtt):

> magn_db:=5;

magn_db := 5

A mágnesek póluserőssége (pozitív számként adjuk meg ezeket, az első az inga):

> magn_p:=
[6.6E-5,

3.6E-5,

3.6E-5,

3.6E-5,

3.6E-5];

# Weber = Wb = V*s

magn_p := [0.66e-4, 0.36e-4, 0.36e-4, 0.36e-4, 0.36e-4]

A mágnesek hossza (a pozitív hossz jelenti, hogy fölfelé áll a pozitív pólus, az első az inga):

> magn_l:=
[+0.06,

-0.06,

-0.06,

-0.06,

-0.06];

# meter = m

magn_l := [0.6e-1, -0.6e-1, -0.6e-1, -0.6e-1, -0.6e-1]

A mágnesek középpontjainak koordinátái (az első az inga):

> magn_choords:=
[Vector([   0.0 ,   0.0 ,    0.012]),

Vector([   0.15,   0.00,   -0.050]),

Vector([   0.05,  -0.00,   -0.050]),

Vector([  -0.05,   0.00,   -0.050]),

Vector([  -0.15,  -0.00,   -0.050])];

# meter = m

magn_choords := [Vector[column]([[0.], [0.], [0.12e-1]], [

Mágnességtani számítások

Az asztal fölötti adott pontban érvényes térerősséget B eljárás számolja ki:

a for ciklus összegzi az asztalhoz rögzített mágnesek hatását B0 eredővektorba,

r1 az i. mágnes pozitív pólusának, r2 a negatív pólusának helyvektora.

> B:=proc(r::Vector)
local B0,B1,B2,r1,r2,i;

B0:=Vector(3,0);

for i from 2 to magn_db do

 r1:=r-magn_choords[i];

 r1[3]:=r1[3]-5/12*magn_l[i];

 B1:=(  magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r1,r1))^3) )  )*r1;

 r2:=r-magn_choords[i];

 r2[3]:=r2[3]+5/12*magn_l[i];

 B2:=-( magn_p[i]  /  ( 4*Pi* (sqrt(DotProduct(r2,r2))^3) )  )*r2;

 B0:=B0+B1+B2;

od;

return B0;

end proc:

Grafikus megjelenítés kiegészítései

A grafikonok színezéséhez itt egyedi függvényt lehet definiálni.

> szin:=proc(a)
return a;

end proc;

szin := proc (a) return a end proc

A feladat fizikai megfogalmazása

A mozgás differenciálegyenletei:

diff(x(t), t, t)*m+diff(x(t), t)*k+x(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[1]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))*magn_p[1])/mu

és

diff(y(t), t, t)*m+diff(y(t), t)*k+y(t)*m*g/l = (B(Vector([x(t), y(t), z+5*magn_l[1]/12]))[2]*magn_p[1]-B(Vector([x(t), y(t), z-5*magn_l[1]/12]))[2]*magn_p[1])/mu

ahol m az inga tömege, l az inga hossza, g a gravivációs gyorsulási együttható, k a közegellenállás.

A megoldások ábrázolása

> m:=0.2; k:=.01; l:=5; g:=9.81; z:=magn_choords[1][3];
me1:=

 diff(x(t),t,t)*m+diff(x(t),t)*k+x(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[1]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[1]*magn_p[1] )/mu:

me2:=

 diff(y(t),t,t)*m+diff(y(t),t)*k+y(t)*m*g/l=

  (B(Vector([x(t),y(t),z+5/12*magn_l[1]]))[2]*magn_p[1]-

   B(Vector([x(t),y(t),z-5/12*magn_l[1]]))[2]*magn_p[1] )/mu:

m := .2

k := 0.1e-1

l := 5

g := 9.81

z := 0.12e-1

> x0:=0.2: y0:=0.0: vx:=-0.1: vy:=0.01: tm:=30: xm:=0.3: ym:=0.3:
a1:=DEplot3d([me1,me2],

      [x(t),y(t)],t=0..tm,x=-xm..xm,y=-ym..ym,

      [[x(0)=x0,y(0)=y0,(D)(x)(0)=vx,(D)(y)(0)=vy]],

      stepsize=0.01,linecolor=szin(t),thickness=2):

list_polys:=[seq([[0 ,magn_choords[i][1],magn_choords[i][2] ],

                 [tm,magn_choords[i][1],magn_choords[i][2] ]],

                i=2..magn_db)]:

a4:=plots[polygonplot3d](list_polys, thickness=3):

plots[display]({a1,a4},orientation=[90,180]);

plots[display]({a1,a4},orientation=[-90,90]);

plots[display]({a1,a4},orientation=[1,88]);

[Plot]

[Plot]

[Plot]

>

Fontos, hogy legyen egy minimális y irányú kezdősebesség, különben az első mágnesen túl rezegne az inga és nem hagyná el az x tengelyt - a maple szerint. A valóságban viszont mingid van akkora egyenetlenség valamiben, ami eldönti merről kerüli meg az első ingát.

A nullánál nagyobb közegellenállás is fontos lehet bizonyos numerikus megoldások esetében, különben lehet, hogy az inga energiája összességében nőne (ami elképzelhetetlen) és "elszállna". Ez a számítási pontatlanságokból adódik - a maple fejlett numerikus technikáinál viszont ilyen probléma nem lép fel.