MODELOWANIE I SYMULACJA

Wyprowadzanie równań różniczkowych ruchu prostych układów mechanicznych

Zadanie 1.

Wyprowadzić równanie różniczkowe opisujące ruch oscylatora harmonicznego wiedząc, że energia
kinetyczna i potencjalna określone są wzorami

2

2

1

1

,

2

2

E

mv

U

kx





gdzie m i k oznaczają odpowiednio masę oscylatora i współczynnik sprężystości.
Rozwiązać otrzymane równanie dla dowolnych warunków początkowych.

Zadanie 2.

Wyprowadzić równanie różniczkowe opisujące ruch wahadła matematycznego wiedząc, że energia
kinetyczna i potencjalna określone są wzorami

2

2

1

,

cos( )

2

E

ml

U

mgl





 



gdzie m oznacza masę, l – długość wahadła, g – przyspieszenie ziemskie.
Rozwiązać otrzymane równanie dla dowolnych warunków początkowych.

Zadanie 3.

Wyprowadzić równania różniczkowe opisujące ruch wahadła matematycznego zawieszonego na
elastycznej lince wiedząc, że energia kinetyczna i potencjalna określone są wzorami









2

2

2

2

2

1

2

1

2

x

y

E

m v

v

U

k

x

y

l

mgy













gdzie m oznacza masę, l – długość wahadła, g – przyspieszenie ziemskie, k – współczynnik
sprężystości.

Zadanie 4.

Wyprowadzić równania różniczkowe opisujące ruch wahadła matematycznego podwieszonego na
oscylatorze wiedząc, że energia kinetyczna i potencjalna określone są wzorami





2

2

2

2

1 1

2

1

1 2

2

2

2

1

2

2

1

1

2

cos(

)

2

2

1

cos(

)

2

E

m v

m v

v v l

x

l v

U

kx

m gl

x













gdzie m

1

oznacza masę oscylatora, m

2

– masę wahadła, l – długość wahadła, g – przyspieszenie

ziemskie, k – współczynnik sprężystości sprężyny.

Wyznaczyć rozwiązanie układu równania za pomocą komendy dsolve z opcją numeric
przyjmując następujące dane: m

1

= 1, m

2

= 0.1, g = 9.81, l = 0.1, k = 100 i warunki początkowe:

1

1

2

2

(0)

0.5 , (0)

0,

(0)

,

(0)

0.

3

x

l v

x

v













Dokonać animacji ruchu wahadła korzystając z komend

> X:=t->'rhs(roz(t)[2])';Y:=t->'rhs(roz(t)[4])';
> a1:=animate(pointplot,[[X(t),0],symbol=box,symbolsize=30],
t=0..5,frames=200,color=blue,color=red):
> a2:=animate(pointplot,[[X(t)+l*sin(Y(t)),-l*cos(Y(t))],
symbol=circle,symbolsize=30],t=0..5,frames=200,color=blue,
color=red):
> a3:=animate(plot,[[[X(t),0],
[X(t)+l*sin(Y(t)),-l*cos(Y(t))]]],t=0..5,axes=none,
scaling=constrained,thickness=2,frames=200):
> display(a1,a2,a3);

>

restart:

>

with(plots):

>

lagrange := proc (n, q, r, L) local i, uzm_q, uzm_r, rel_r_q, Lq,

Lr, Lrt; global row; uzm_q := seq(q[i] = q[i](t), i = 1 .. n); uzm_r
:= seq(r[i] = r[i](t), i = 1 .. n); for i to n do Lq[i] := subs([uzm_q,
uzm_r], diff(L, q[i])); Lr[i] := subs([uzm_q, uzm_r], diff(L, r[i]))
end do; for i to n do Lrt[i] := diff(Lr[i], t) end do; rel_r_q :=
seq(r[i](t) = diff(q[i](t), t), i = 1 .. n); for i to n do row[i] :=
subs(rel_r_q, Lrt[i]-Lq[i] = 0) end do; seq(row[i], i = 1 .. n) end
proc:

>
>

#Zad1

>

E := (1/2)*m*v[1]^2:

>

U := (1/2)*k*x^2:

>

L := E-U:

>

lagrange(1, x, v, L):

>
>

#Zad2

>

restart:

>

E := (1/2)*m*l^2*omega[1]^2:

>

U := -m*g*l*cos(phi[1]):

>

L := E-U:

>

lagrange(1, phi, omega, L):

>
>

#Zad3

>

restart:

>

E := (1/2)*m*(v[1]^2+v[2]^2):

>

U := (1/2)*k*(sqrt(x[1]^2+x[2]^2)-l)^2+m*g*x[2]:

>

L := E-U:

>

lagrange(2, v, x, L):

>
>

#Zad4

>

restart:

>

E :=

(1/2)*m[1]*v[1]^2+(1/2)*m[2]*(v[1]^2+2*v[1]*v[2]*l*cos(x[2])+l^2*
v[2]^2):

>

U := (1/2)*k*x[1]^2-m[2]*g*l*cos(x[2]):

>

m[1] := 1: m[2] := .1: g := 9.81: l := .1: k := 100:

>

wp := x[1](0) = .5*l, (D(x[1]))(0) = 0, x[2](0) = (1/3)*Pi,

(D(x[2]))(0) = 0:

>

L := E-U:

>

lagrange(2, x, v, L):

>

dsolve({wp, row[1], row[2]}, {x[1](t), x[2](t)}, numeric):

MODELOWANIE I SYMULACJA

Modelowanie liniowego układu dyskretnego

Wyznaczyć rozwiązanie analityczne opisujące ruch liniowego układu dyskretnego przedstawionego
na rysunku

gdzie m oznacza masę poszczególnych ciał, k – sztywność każdej sprężyny, a x

i

, i = 1, 2, ...,5 –

przemieszczenia środków poszczególnych mas.

Kolejność postępowania:

1. Wyznaczyć energię kinetyczną i potencjalną układu wykorzystując zależności:

5

4

2

2

1

1

1

1

1

,

(

)

2

2

i

i

i

i

i

E

m

v

U

k

x

x

















2. Wyznaczyć macierz sztywności K i bezwładności M korzystając ze wzorów:

2

2

,

ij

ij

i

j

i

j

U

E

k

m

x x

v v









 

 

Uwaga: należy zadeklarować najpierw obie macierze, by następnie wyznaczyć ich elementy
w podwójne pętli

3. Wyznaczyć macierz A związaną z macierzami K i M wzorem





1

A

M K

4. Przyjąć dane liczbowe: m = 1, k = 100 i wyznaczyć wartości i wektory własne macierzy A

korzystając z komendy Eigenvectors w postaci

> alpha,V:=Eigenvectors(A):

5. Korzystając z procedury sortowanie posortować wartości własne zmieniając

równocześnie w odpowiedni sposób pozycje kolumn w macierzy V, zawierającej wektory
własne macierzy A

> alpha,V:=sortowanie(alpha,V);

6. Dokonać normalizacji wektorów własnych względem macierzy bezwładności korzystając ze

wzoru (u

i

– kolumny macierzy V)





,

1, 2,...,

i

i

T
i

i

i

n





u

w

u Mu

7. Wyznaczyć częstości drgań swobodnych korzystając z komendy map ze wzoru

,

1, 2,...,

i

i

i

n

  



gdzie



i

oznaczają posegregowane wartości własne.

8. Zadać warunki początkowe w formie wektora przemieszczeń początkowych x0

zawierającego zerowe elementy i wektora prędkości początkowych v0 zawierającego
niezerową pierwszą współrzędną równą 10

> x0:=Vector(5,[0,0,0,0,0]); v0:=Vector(5,[10,0,0,0,0]);

9. Wyznaczyć analityczne rozwiązanie opisujące ruch poszczególnych mas korzystając ze

wzoru

5

1

1

2

0

. .( 0

0 )

. .

0 cos (

)

sin (

)

T

T
i

i

i

i

i

i

t

t

t



















 





 































v

x

w M x

v

w

w M x

w

10. Sporządzić wykres przemieszeń i prędkości poszczególnych mas korzystając z poniższych

komend:

> plot([seq(x[i],i=1..n)],t=0..3);
> v:=map(diff,x,t):
> plot([seq(v[i],i=1..n)],t=0..3);

11. Dokonać animacji ruchu poszczególnych mas korzystając z poniższych komend:

> seq_pkt:=seq([(j-1)*3+'x'[j],0],j=1..n);
> animate(pointplot,[[seq_pkt],symbol=box,symbolsize=50],
t=0..3,frames=100,color=red,axes=none);

12. Sporządzić wykres czasowy energii mechanicznej (E + U)

> plot(E+U,t=0..5,0..100);

13. Która ze sprężyn będzie najbardziej ściśnięta/rozciągnięta w trakcie pierwszych trzech

sekund?

>

restart:

>

with(LinearAlgebra):

with(plots):

>

#sortowanie proc

>
>

#Zad1

>

E:=1/2*m*add(v[i]^2,i=1..5):

>

U:=1/2*k*add((x[i+1]-x[i])^2,i=1..4):

>
>

#Zad2

>

n:=5:

>

K:=Matrix(n):

>

M:=Matrix(n):

>

for i to n do

for j to n do
K[i,j]:=diff(U,x[i],x[j]);
M[i,j]:=diff(E,v[i],v[j]);
end do;
end do;

>

K:

>

M:

>
>

#Zad3

>

A:=1/M.K:

>
>

#Zad4

>

m:=1: k:=100:

>

alpha,V:=evalf(Eigenvectors(A)):

>
>

#Zad5

>

alpha,V:=sortowanie(alpha,V):

>
>

#Zad6

>

for i from 1 to n do

w[i]:=V[1..n,i]/(sqrt((V[1..n,i]^%T).M.V[1..n,i])):
end do:

>
>

#Zad7

>

for i from 1 to n do

omega[i]:=sqrt(alpha[i])
end do:

>
>

#Zad8

>

x0:=Vector(5,[0,0,0,0,0]):

>

v0:=Vector(5,[10,0,0,0,0]):

>
>

#Zad9

>

x:=(Transpose(w[1]).M.(x0+v0*t))*w[1]+add((Transpose(w[1]).M.(x0*cos(om
ega[i]*t)+v0*sin(omega[i]*t)/omega[i]))*w[i], i=2..5):