Difference between revisions of "Махало"

From Ilianko
 
(39 intermediate revisions by the same user not shown)
Line 13: Line 13:
 
<math> s = \ell \theta </math>
 
<math> s = \ell \theta </math>
  
За да изразим скоростта ползваме първата производна (изменението на разстоянието за единица време). Изменя се само ъгълът, а дължината на нишката е константа и се запазва:
+
За да изразим скоростта ползваме първата [[производна]] (изменението на разстоянието за единица време). Изменя се само ъгълът, а дължината на нишката е константа и се запазва:
  
 
<math> v = {ds\over dt} = {{d\ell\theta}\over dt}=  \ell {d\theta\over dt}</math>
 
<math> v = {ds\over dt} = {{d\ell\theta}\over dt}=  \ell {d\theta\over dt}</math>
Line 20: Line 20:
  
 
<math> a = {d^2s\over dt^2} = \ell{d^2\theta\over dt^2} </math>  (1)
 
<math> a = {d^2s\over dt^2} = \ell{d^2\theta\over dt^2} </math>  (1)
 +
 
== Сили действащи на махалото ==
 
== Сили действащи на махалото ==
  
 
Под внимание се взема само силата действаща по посока на движението. Допирателната по окръжността на движение (тангенциалната сила ( танго -> допир)). Перпендикулярната сила (нормалната ) се неутрализира.
 
Под внимание се взема само силата действаща по посока на движението. Допирателната по окръжността на движение (тангенциалната сила ( танго -> допир)). Перпендикулярната сила (нормалната ) се неутрализира.
  
От II закон на нютон <math> F = ma\,</math>
+
От [[II закон на Нютон]] <math> F = ma\,</math>
  
 
За махалото <math> F = gm \sin\theta = -ma\, => a = -g\sin\theta\, </math>  (2)
 
За махалото <math> F = gm \sin\theta = -ma\, => a = -g\sin\theta\, </math>  (2)
Line 37: Line 38:
 
<math>\theta < 20^{\circ}</math>. В този случай диференциалното уравнение ще се преобразува в линейно от втори ред.
 
<math>\theta < 20^{\circ}</math>. В този случай диференциалното уравнение ще се преобразува в линейно от втори ред.
  
=== Аналаитично решение ===
+
=== Аналитично решение ===
  
 
<math>\ell{d^2\theta\over dt^2} + g\sin\theta = 0  <=> \ell{d^2\theta\over dt^2} + g\theta = 0, \theta < 20^{\circ}</math>
 
<math>\ell{d^2\theta\over dt^2} + g\sin\theta = 0  <=> \ell{d^2\theta\over dt^2} + g\theta = 0, \theta < 20^{\circ}</math>
  
Полагаме <math>\ {\theta} = {y} </math>
+
Полагаме <math>\ {\theta} = {y} </math>, от тук <math> \ell{y''} + {g}{y} = 0 </math> - хомогенно диференциално уравнение от 2-ри ред
 +
 
 +
Уравнението е от вида <math>\ ay''+by'+cy = 0</math>, затова лесно могат да се намерят решения във вида
 +
<math>\ y = c_{1} e^{r_{1} x} + c_2 e^{r_{2} x} </math>, където за ''r'' се решава уравнението <math> \ ar^{2}+br+c=0</math>, a именно:
 +
 
 +
<math> \ell r^2 + g = 0 </math>
 +
 
 +
<math> r = \pm j \sqrt{\frac{g}{l}} </math>
 +
 
 +
<math>\ y(t) =  c_{1} e^{-j \sqrt{\frac{g}{l}} t} + c_2 e^{j \sqrt{\frac{g}{l}} t}</math>
 +
 
 +
<math>\ e^{jx} = cos(x) + j sin(x) </math> от тук =>
 +
 
 +
<math>\ \theta(t) = C_{1}cos(\sqrt{\frac{g}{l}} t) + C_{2}sin(\sqrt{\frac{g}{l}} t) </math>
 +
 
 +
За получаване на конкретно решение ще зададем начални условия. Примерно в началното положение <math> \ t=0 </math>, махалото е отклонено на <math>\ \theta_{0} </math> градуса и пуснато свободно, т.е. с нулева начална скорост.
 +
 
 +
<math>
 +
\begin{cases}
 +
\theta(0) = \theta_{0} = C_{1}cos(\sqrt{\frac{g}{l}} 0) + C_{2}sin(\sqrt{\frac{g}{l}} 0) \\
 +
\theta(0)' = 0 = -C_{1}\sqrt{\frac{g}{l}}sin(\sqrt{\frac{g}{l}} 0) + C_{2}\sqrt{\frac{g}{l}}cos(\sqrt{\frac{g}{l}} 0)
 +
\end{cases} </math>
 +
=>
 +
<math>\ C_{1} = \theta_{0} , C_{2} = 0
 +
 
 +
</math> 
 +
 
 +
За решение на задача при малки ъгли се получава
 +
 
 +
<math> \theta(t) = \theta_{0} cos(\sqrt{\frac{g}{l}} t)  </math>
 +
 
 +
<math> \sqrt{\frac{g}{l}} = \sqrt{\frac{m}{s^{2}}*\frac{1}{m}} = \frac{1}{s} = \omega * \frac{rad}{s} </math> e ъгловата скорост с която ще се движи махалото.
 +
 
 +
От тук <math>\ \theta(t) = \theta_{0} cos(\omega t)  </math>
 +
 
 +
=== Визуализация на решението с Matlab ===
 +
 
 +
Имаме следните начални условия:
 +
*<math> \ l = 5  </math>
 +
*<math> \theta_{0} = 17^{\circ} </math>
 +
 
 +
<code><pre>
 +
l = 5; % m
 +
g = 9.8; % m/s^2
 +
theta_0 = 17; % degree
 +
theta_0r = 17*pi/180 % ъгъл в радиани
 +
omega = sqrt(g/l); % ъгловата скорост
 +
t = 0:0.1:20; % интуитивно си избираме подходящ интервал за изследване
 +
theta = theta_0r*cos(omega.*t); % Фундаментална формула, с който започва анализа на трептенията
 +
 
 +
% За тези с по-малко въображение и тези, които искат да представят резултатите на други
 +
plot(t,theta);
 +
xlabel('t [s]'), ylabel('theta [rad]');
 +
axis([ 0 20 -0.4 0.4]);
 +
title('Matematichesko mahalo');
 +
% Само трябва да добавим съпротивление на въздуха; вятър;
 +
% маса на нишката, която държи махалото; трептене на опoрната точка и да решим за всякакви начални условия.
 +
</pre></code>
 +
 
 +
[[Image:mahAn.png|frame|none|center|Нито повече, нито по-малко от една синосоида]]
 +
 
 +
=== Аналитично решение с Matlab ===
 +
Общо решение:
 +
<code><pre>
 +
syms g l;
 +
y = dsolve('D2y = -(g/l)*y', 't')
 +
 
 +
%решение
 +
% yC15*exp((t*(-g*l)^(1/2))/l) + C16/exp((t*(-g*l)^(1/2))/l)
 +
 
 +
</pre></code>
 +
 
 +
Частно решение с начални условия
 +
<code><pre>
 +
% l=5, g=9.8, 0.3 = 17*pi*180
 +
y = dsolve('D2y = -(9.8/5)*y','Dy(0)=0, y(0) =0.3','t')
 +
 
 +
%решение
 +
% y = (3*cos((7*t)/5))/10
 +
</pre></code>
 +
 
 +
=== Числено решаване с Matlab ===
 +
 
 +
Решението се състои в две части.
 +
Първо се създава файл - функция описваща уравнението
 +
Второ изпълнява се функцията ode45();
 +
 
 +
''' Функция mah1 '''
 +
<code><pre>
 +
function  ds = mah1(t, yy, l , g)
 +
% t е променливата, по която диференцираме, в случая не се ползва
 +
% защото не е част от уравнението
 +
% y функцията, която търсим
 +
% dy  e  нейната първа производна
 +
% d^2y/dt^2 + g/l*sin(y) = 0 - matematichesko mahalo
 +
 
 +
y = yy(1);  %  Тези два реда са за онагледяване,
 +
dy = yy(2); %  стойностите за уу може дирктно да се
 +
% запишат в долните изрази.
 +
 
 +
ds(1,1) = dy;  % <=> ( y(t) )' = dy(t)
 +
ds(2,1) = -g/l*sin(y); % <=> ( dy(t) )' = -g/l*sin(y)
 +
% ds са конкретни числови стойности
 +
end
 +
</pre></code>
 +
 
 +
''' Скрипт за извикване на ode45()
 +
<code><pre>
 +
g=9.8;
 +
l=5;
 +
x0 = [17*pi/180; 0]; % начални условия
 +
t = 0:0.01:20;
 +
options = [] ; % за сега не ни интересуват
 +
[t,x] = ode45(@mah1,t,x0,options,l,g);
 +
 
 +
plot(t, x(:,1)); % make the plot
 +
title('Reshenie mat. mahalo s ode45()');
 +
</pre></code>
 +
 
 +
== [[Числено интегриране]] ==
 +
При численото решаване, знаем началната стойност на функцията и [[производна|начина по който тя се изменя]]. С тези входни данни, точка по точка се изчислява всяка следваща стойност на функцията.
 +
 
 +
Най-лесният алгоритъм за тази цел е [http://en.wikipedia.org/wiki/Euler_method методът на Ойлер].
 +
Ето какво образно показано изпълнява ode45()
 +
 
 +
<code><pre>
 +
% Метод на Ойлер за интегриране
 +
g= 9.8; l = 5; % константи
 +
T = 20; % максимална стойност, до която ще интегрираме
 +
dt = 0.0001;  % колкото по-малко, толкова по-точно и по-бавно
 +
t = 0:dt:T;  % интервала, който ще разгледаме
 +
points = T/dt;
 +
 
 +
%масив, където ще се съхранят всички точки
 +
ydot = zeros(1,points+1);
 +
y = zeros(1,points+1);
 +
 
 +
% Начални условия
 +
y(1) = 17*pi/180;
 +
ydot(1) = 0;
 +
 
 +
for i = 1:T/dt
 +
    y(i+1) = y(i) + ydot(i)*dt;
 +
    ydot(i+1) = ydot(i) - g/l*sin(y(i))*dt;
 +
end
 +
plot(t, y)
 +
</pre></code>
 +
 
 +
Друг вариант
 +
<code><pre>
 +
g= 9.8; l = 5; % константи
 +
order = 2;
 +
T = 20;
 +
dt = 0.0001
 +
points = T/dt;
 +
%y = zeros(1,points);
 +
%t = zeros(1,points);
 +
 
 +
x(1) = 17*pi/180;
 +
x(2) = 0;
 +
 
 +
for i = 1:points
 +
   
 +
    fxt(1) = x(2);
 +
    fxt(2) = - g/l*sin(x(1));
 +
   
 +
    for count = 1:order
 +
        x(count) = x(count) + dt*fxt(count);
 +
    end
 +
  y(i+1) = x(1);
 +
  t(i+1) = t(i) + 1;
 +
end
 +
plot(t, y)
 +
</pre></code>
 +
 
 +
И вариант на C
 +
<code><pre>
 +
#include <stdio.h>
 +
#include <math.h>
 +
 
 +
/* Variables defined here are global.  */
 +
float x[2],xdot[2];    /* Define states and state derivatives  */
 +
float t,dt;            /* Time and integration interval global */
 +
int nstate;            /* Number of states and state equations */
 +
float g = 9.8, l = 5;
 +
 
 +
/* State Equations */
 +
int state_eqns(void)
 +
{
 +
    xdot[0] = x[1];
 +
    xdot[1] = -g/l*sin(x[0]) ;
 +
    return 0;
 +
}
 +
 
 +
/* Euler Integration Algorithm */
 +
void euler(void)
 +
{
 +
int count = 0;
 +
    state_eqns();
 +
    while (count<nstate)
 +
    {
 +
x[count] = x[count] + dt*xdot[count];
 +
        count = count+1;
 +
    }
 +
    t = t + dt;
 +
    return;
 +
}
 +
 
 +
int main(void)
 +
{
 +
int    index,iterations;
 +
    float  Tmax;
 +
 
 +
    /* Data Input Section *******************************************/
 +
    printf(" Input data for the simulation. \n");
 +
    printf("\n");
 +
   
 +
    printf("Input the total time for your simulation.\n");
 +
    scanf("%f",&Tmax);
 +
   
 +
    printf("Input the sample time for your simulation.\n");
 +
    scanf(" %f",&dt);
 +
 
 +
    /* Initialization section ***************************************/
 +
    t=0.;          /* Initialize time & states*/
 +
    index=0;
 +
    nstate=2;  /* First order system - 1 state */
 +
         
 +
    x[0] = 0.2967; /* 17*pi/180 */
 +
    xdot[0] = 0;
 +
    x[1] = 0.;
 +
    xdot[1] = 0.;
 +
   
 +
    index = 0;
 +
    iterations = Tmax/dt;
 +
 
 +
    /* Do simulation calculations and print results ************************/
 +
    while (index<=iterations)
 +
    {
 +
euler();
 +
        printf("t(%i)=%12.4f; m(%i)=%12.4f; \n",index,t,index,x[0]);
 +
        index = index + 1;
 +
    }
 +
 
 +
    return 0;
 +
}
 +
</pre></code>
 +
 
 +
'''Вариант с PHP'''
 +
<code><pre>
 +
 
 +
<?
 +
// Метод на Ойлер за интегриране
 +
 
 +
$g= 9.8; $l = 5; //'константи';
 +
$T = 20; //максимална стойност, до която ще интегрираме
 +
$dt = 0.001;  // колкото по-малко, толкова по-точно и по-бавно
 +
 
 +
 
 +
//%масив, където ще се съхранят всички точки
 +
$ydot = 0;
 +
$y = array(0,0);
  
<math> \ell{y''} + {g}{y} = 0 </math>
+
//% Начални условия
 +
$y[1] = 17*pi()/180;
 +
$ydot = 0;
  
 +
for ($i = 1; $i <= $T/$dt; $i++)
 +
{
 +
  array_push($y,($y[$i] + $ydot*$dt));
 +
  $ydot = $ydot - $g/$l*sin($y[$i])*$dt;
 +
}
 +
 
 +
print_r($y);
  
 +
?>
 +
</pre></code>
  
 +
== Връзки ==
 +
[http://dw.georgievi.net/ivan/ Бифилярно махало, Пружинно махало, Физично махало]
  
  

Latest revision as of 10:10, 5 January 2013

Движение на махало под въздействие на гравитацията

Цел: Създаване на връзки между диференциране и диференциални уравнения.

Постановка: Математическо махало. (Съпротивлението на въздуха се пренебрегва). Да се намери зависимостта на махалото спрямо времето, т.е. функцията на движение на махалато.

Анализ: Във всеки момент махалото се движи с различна скорост и ускорение, зависещи от ъгъла на отклонение на махалото.

Производни, разстояние, скорост и ускорение.

Движението на махалото се извършва по окръжност с радиус равен на дължината на нишката, на която е окачено махалото. От тук, следва че:

За да изразим скоростта ползваме първата производна (изменението на разстоянието за единица време). Изменя се само ъгълът, а дължината на нишката е константа и се запазва:

По да изразим ускорението, използваме втората производна

(1)

Сили действащи на махалото

Под внимание се взема само силата действаща по посока на движението. Допирателната по окръжността на движение (тангенциалната сила ( танго -> допир)). Перпендикулярната сила (нормалната ) се неутрализира.

От II закон на Нютон

За махалото (2)

Диференциално уравнение на махало

От (1) и (2)

Решаване на диференциалното уравнение

Има два подхода при решаване на диференциални уравнения, аналитичен и числен. При конкретното уравнение са възможни и двата подхода, но аналитичният метод може да се използва за . В този случай диференциалното уравнение ще се преобразува в линейно от втори ред.

Аналитично решение

Полагаме , от тук - хомогенно диференциално уравнение от 2-ри ред

Уравнението е от вида , затова лесно могат да се намерят решения във вида , където за r се решава уравнението , a именно:

от тук =>

За получаване на конкретно решение ще зададем начални условия. Примерно в началното положение , махалото е отклонено на градуса и пуснато свободно, т.е. с нулева начална скорост.

=>

За решение на задача при малки ъгли се получава

e ъгловата скорост с която ще се движи махалото.

От тук

Визуализация на решението с Matlab

Имаме следните начални условия:

l = 5; % m
g = 9.8; % m/s^2
theta_0 = 17; % degree
theta_0r = 17*pi/180 % ъгъл в радиани
omega = sqrt(g/l); % ъгловата скорост
t = 0:0.1:20; % интуитивно си избираме подходящ интервал за изследване
theta = theta_0r*cos(omega.*t); % Фундаментална формула, с който започва анализа на трептенията

% За тези с по-малко въображение и тези, които искат да представят резултатите на други
plot(t,theta); 
xlabel('t [s]'), ylabel('theta [rad]');
axis([ 0 20 -0.4 0.4]);
title('Matematichesko mahalo');
% Само трябва да добавим съпротивление на въздуха; вятър;
% маса на нишката, която държи махалото; трептене на опoрната точка и да решим за всякакви начални условия. 
Нито повече, нито по-малко от една синосоида

Аналитично решение с Matlab

Общо решение:

syms g l;
y = dsolve('D2y = -(g/l)*y', 't')

%решение
% yC15*exp((t*(-g*l)^(1/2))/l) + C16/exp((t*(-g*l)^(1/2))/l)

Частно решение с начални условия

% l=5, g=9.8, 0.3 = 17*pi*180
y = dsolve('D2y = -(9.8/5)*y','Dy(0)=0, y(0) =0.3','t') 

%решение
% y = (3*cos((7*t)/5))/10

Числено решаване с Matlab

Решението се състои в две части. Първо се създава файл - функция описваща уравнението Второ изпълнява се функцията ode45();

Функция mah1

function  ds = mah1(t, yy, l , g)
% t е променливата, по която диференцираме, в случая не се ползва
% защото не е част от уравнението
% y функцията, която търсим
% dy  e  нейната първа производна
% d^2y/dt^2 + g/l*sin(y) = 0 - matematichesko mahalo

y = yy(1);  %  Тези два реда са за онагледяване,
dy = yy(2); %  стойностите за уу може дирктно да се
% запишат в долните изрази.

ds(1,1) = dy;  % <=> ( y(t) )' = dy(t)
ds(2,1) = -g/l*sin(y); % <=> ( dy(t) )' = -g/l*sin(y)
% ds са конкретни числови стойности
end

Скрипт за извикване на ode45()

g=9.8;
l=5;
x0 = [17*pi/180; 0]; % начални условия
t = 0:0.01:20;
options = [] ; % за сега не ни интересуват
[t,x] = ode45(@mah1,t,x0,options,l,g);

plot(t, x(:,1)); % make the plot
title('Reshenie mat. mahalo s ode45()');

Числено интегриране

При численото решаване, знаем началната стойност на функцията и начина по който тя се изменя. С тези входни данни, точка по точка се изчислява всяка следваща стойност на функцията.

Най-лесният алгоритъм за тази цел е методът на Ойлер. Ето какво образно показано изпълнява ode45()

% Метод на Ойлер за интегриране
g= 9.8; l = 5; % константи
T = 20; % максимална стойност, до която ще интегрираме
dt = 0.0001;   % колкото по-малко, толкова по-точно и по-бавно
t = 0:dt:T;  % интервала, който ще разгледаме
points = T/dt;

%масив, където ще се съхранят всички точки
ydot = zeros(1,points+1);
y = zeros(1,points+1);

% Начални условия
y(1) = 17*pi/180;
ydot(1) = 0;

for i = 1:T/dt
    y(i+1) = y(i) + ydot(i)*dt;
    ydot(i+1) = ydot(i) - g/l*sin(y(i))*dt; 
end
plot(t, y)

Друг вариант

g= 9.8; l = 5; % константи
order = 2;
T = 20;
dt = 0.0001
points = T/dt;
%y = zeros(1,points);
%t = zeros(1,points);

x(1) = 17*pi/180;
x(2) = 0;

for i = 1:points
    
    fxt(1) = x(2);
    fxt(2) = - g/l*sin(x(1)); 
    
    for count = 1:order
        x(count) = x(count) + dt*fxt(count);
    end
   y(i+1) = x(1);
   t(i+1) = t(i) + 1;
end
plot(t, y)

И вариант на C

#include <stdio.h> 
#include <math.h>

/* Variables defined here are global.  */ 
float x[2],xdot[2];     /* Define states and state derivatives  */ 
float t,dt;             /* Time and integration interval global */ 
int nstate;             /* Number of states and state equations */ 
float g = 9.8, l = 5; 

/* State Equations */ 
int state_eqns(void) 
{ 
    xdot[0] = x[1];
    xdot[1] = -g/l*sin(x[0]) ;
    return 0; 
} 

/* Euler Integration Algorithm */ 
void euler(void) 
{ 
	int count = 0; 
    state_eqns(); 
    while (count<nstate) 
    { 
		x[count] = x[count] + dt*xdot[count];
        count = count+1; 
    } 
    t = t + dt;
    return; 
} 

int main(void) 
{ 
	int     index,iterations; 
    float  Tmax; 

    /* Data Input Section *******************************************/ 
    printf(" Input data for the simulation. \n"); 
    printf("\n"); 
    
    printf("Input the total time for your simulation.\n"); 
    scanf("%f",&Tmax); 
    
    printf("Input the sample time for your simulation.\n"); 
    scanf(" %f",&dt); 

    /* Initialization section ***************************************/ 
    t=0.;          /* Initialize time & states*/ 
    index=0; 
    nstate=2;  /* First order system - 1 state */ 
          
    x[0] = 0.2967; /* 17*pi/180 */ 
    xdot[0] = 0;
    x[1] = 0.;
    xdot[1] = 0.;
    
    index = 0; 
    iterations = Tmax/dt; 

    /* Do simulation calculations and print results ************************/ 
    while (index<=iterations) 
    { 
		euler(); 
        printf("t(%i)=%12.4f; m(%i)=%12.4f; \n",index,t,index,x[0]); 
        index = index + 1; 
    } 

    return 0; 
}

Вариант с PHP


<?
// Метод на Ойлер за интегриране

$g= 9.8; $l = 5; //'константи';
$T = 20; //максимална стойност, до която ще интегрираме
$dt = 0.001;   // колкото по-малко, толкова по-точно и по-бавно


//%масив, където ще се съхранят всички точки
$ydot = 0;
$y = array(0,0);

//% Начални условия
$y[1] = 17*pi()/180;
$ydot = 0;

for ($i = 1; $i <= $T/$dt; $i++)
{
   array_push($y,($y[$i] + $ydot*$dt));
   $ydot = $ydot - $g/$l*sin($y[$i])*$dt;
}
  
print_r($y);

?>

Връзки

Бифилярно махало, Пружинно махало, Физично махало