Задание2

задание2.1. С помощью инструментария символьной математики (Symbolic Toolbox) построить аналитическое решение, описывающее свободные колебания системы при наличии вязкого трения.

>> y=dsolve('D2y+2*eps*D1y+theta^2*y=0','y(0)=y0','Dy(0)=v0','t')

y =

-1/2*(eps^2-theta^2)^(1/2)*((eps^2-theta^2)^(1/2)*y0+v0+eps*y0)/(-eps^2+theta^2)*exp((-eps+(eps^2-theta^2)^(1/2))*t)+1/2*(eps^2-theta^2)^(1/2)/(-eps^2+theta^2)*(v0+eps*y0-(eps^2-theta^2)^(1/2)*y0)*exp((-eps-(eps^2-theta^2)^(1/2))*t)

>> subs(y,{'y0','v0','theta','eps'},{1,0,3,0.1})

ans =

-1/1798*(-899)^(1/2)*100^(1/2)*(1/10+1/100*(-899)^(1/2)*100^(1/2))*exp((-1/10+1/100*(-899)^(1/2)*100^(1/2))*t)+1/1798*(-899)^(1/2)*100^(1/2)*(1/10-1/100*(-899)^(1/2)*100^(1/2))*exp((-1/10-1/100*(-899)^(1/2)*100^(1/2))*t)

>> ezplot((ans),[0 20 ])

задание 2.2. Численно проинтегрировать уравнение свободных колебаний механической системы при наличии вязкого трения

создаем в редакторе

function task2_2

% формирование вектора начальных условий

Y0 = [1; 0];

% вызов солвера от функции oscil, начального и конечного

% момента времени и вектора начальных условий

theta=3;

eps=0.1;

b=0.0;

[T, Y] = ode15s(@oscil, [0 20], Y0,[],theta,eps,b);

% вывод графика решения исходного дифференциального уравнения

plot(T, Y(:,1)) %plot(T, 3*Y(:,1)+4*Y(:,2))

% подфункция вычисления правых частей уравнений

function F = oscil(t, y,theta,eps,b)

F = [y(2); -b*theta^2*sign(y(2))-2*eps*y(2)-theta^2*y(1)];

и выполняем ее в Command Window, изображая график в том же окне

>> hold on

>> task2_2

графики практически неразличимы!

Задание 2.3 Численно проинтегрировать уравнение свободных колебаний механической системы при наличии сухого трения и комбинации вязкого и сухого трения

модифицируем task2_2 на task2_3, полагая

для СУХОГО трения eps=0.0;b=0.01;

>> task2_3

убывание размахов по линейному закону! верно!!

при eps=0.01;b=0.01 имеет место комбинация трений:

>> task2_3

Warning: Failure at t=4.292833e+001. Unable to meet integration tolerances

without reducing the step size below the smallest value allowed

(1.136868e-013) at time t.

> In ode15s at 751

In task2_3 at 9

возникают проблемы с численным интегрированием при больших временах…

а все дело в том, что останавливается система - в зоне застоя

в самом деле..при eps=0.01;b=0.1;

не смогли построить решение даже на отрезке [0 20]:

т.к. было [T, Y] = ode15s(@oscil, [0 20], Y0,[],theta,eps,b);

Задание 3

Задание 3.1. Построить и визуализировать средствами Symbolic Toolbox аналитическое решение задачи о вынужденных колебаниях при гармоническом воздействии с учетом вязкого сопротивления.

например, такой скрипт

%task3_1

y=dsolve('D2y+2*eps*D1y+theta^2*y=P/M*cos(w*t)','y(0)=y0','Dy(0)=v0','t')

subs(y,{'y0','v0','w','theta','P','M','eps'},{0,0,3,4,5,6,0.1})

ezplot(ans,[0 60 ])

тогда получим

>> task3_1 % может быть стоило поставить; …)))

y =

1/2*exp((-eps+(eps^2-theta^2)^(1/2))*t)*(-y0*M*theta^4+y0*M*theta^2*w^2-eps*v0*M*theta^2+(eps^2-theta^2)^(1/2)*v0*M*theta^2+P*theta^2-2*y0*M*eps^2*w^2-eps*v0*M*w^2-2*eps*w^2*(eps^2-theta^2)^(1/2)*y0*M-(eps^2-theta^2)^(1/2)*v0*M*w^2)*((eps^2-theta^2)^(1/2)*eps-eps^2+theta^2)/M/(w^2-theta^2+2*eps^2-2*(eps^2-theta^2)^(1/2)*eps)/theta^2/(-eps^2+theta^2)+1/2*exp((-eps-(eps^2-theta^2)^(1/2))*t)*(-(eps^2-theta^2)^(1/2)*eps-eps^2+theta^2)*(-eps*v0*M*w^2-eps*v0*M*theta^2+2*eps*w^2*(eps^2-theta^2)^(1/2)*y0*M-2*y0*M*eps^2*w^2+P*theta^2+(eps^2-theta^2)^(1/2)*v0*M*w^2-(eps^2-theta^2)^(1/2)*v0*M*theta^2+y0*M*theta^2*w^2-y0*M*theta^4)/theta^2/(-eps^2+theta^2)/M/(w^2+2*eps^2-theta^2+2*(eps^2-theta^2)^(1/2)*eps)-P*(-cos(w*t)*theta^2+cos(w*t)*w^2-2*sin(w*t)*w*eps)/M/(4*eps^2*w^2+theta^4-2*theta^2*w^2+w^4)

ans =

125/4797*exp((-1/10+1/100*(-1599)^(1/2)*100^(1/2))*t)*(1/1000*(-1599)^(1/2)*100^(1/2)+1599/100)/(-349/50-1/500*(-1599)^(1/2)*100^(1/2))+125/4797*exp((-1/10-1/100*(-1599)^(1/2)*100^(1/2))*t)*(-1/1000*(-1599)^(1/2)*100^(1/2)+1599/100)/(-349/50+1/500*(-1599)^(1/2)*100^(1/2))+875/7404*cos(3*t)+25/2468*sin(3*t)

Задание 3.2. Реализовать решение задачи о колебаниях при воздействии произвольного вида с учетом вязкого сопротивления и при заданных значениях параметров в форме интеграла Дюамеля, вычислив его аналитически и численно

большого смысла в реализации интеграла Дюамеля нет, если есть ОДЕ- солверы..в основном, чтобы потренироваться в интегрировании на МАТЛАБе…

скрипт,реализующий при нулевых НУ интеграл АНАЛИТИЧЕСКИ

% sym_duamel

syms t x w theta eps P M

eta=sqrt(theta^2-eps^2);

f=P/M/eta*cos(w*x)*exp(-eps*(t-x))*sin(eta*(t-x));

I=int(f,x,0,t);

ezplot(char(subs(I,{w, theta, eps, P, M},{3,4,0.1,5,6})),[0 60])

решение неотличимо от предыдущего…

будут проблемы, если интеграл АНАЛИТИЧЕСКИ не вычисляется…поэтому

пытаемся реализовать численно

сначала- функцию, реализующую интеграл с переменным верхним пределом…

function f=duamel(t)

% это основная функция

f=quadl(@fint,0,t,1e-7,0,t);

% а это подфункция

function f=fint(x,t)

P=5;

M=6;

w=3;

eps=0.1;

theta=4;

eta=sqrt(theta^2-eps^2);

f=P/M/eta*cos(w*x).*exp(-eps*(t-x)).*sin(eta*(t-x));

подфункция в этом же файле реализует подынтегральное выражение…

стоят.* все для quadl….

теперь можно запустить

>>fplot(@duamel,[0.01 60])

считает довольно долго и сухое трение не учесть- не стоит использовать интеграл Дюамеля …но результат – тот же..

Задание 3.3 Аналитически с использованием концепции эквивалентного вязкого демпфирования и численно исследовать вынужденные колебания при гармоническом воздействии с учетом сухого трения и комбинации сухого и вязкого трения

Численное решение строится очевидным образом- в подфункцию добавляем Р*соs….

Аналитическое решение с использованием концепции эквивалентного вязкого демпфирования справедливо только для установившегося процесса и сводится просто в к пересчету eps=eps/(1-4*b*theta^2*M/P/pi); и затем к обычному использованию dsolve.

Поэтому реализация задания может быть такой

function task3_3

% формирование вектора начальных условий

Y0 = [1; 2];

% задание параметров и вызов солвера от подфункции oscil

theta=4;

eps=0.05;

b=0.01;

P=5;

M=6;

w=3;

[T, Y] = ode15s(@oscil, [0 80], Y0,[],theta,eps,b,w,P,M);

% вывод графика решения исходного дифференциального уравнения

plot(T, Y(:,1))

% аналитическое решение

eps=eps/(1-4*b*theta^2*M/P/pi);% вычисляем эквивалентный eps

figure

y=dsolve('D2y+2*eps*D1y+theta^2*y=P/M*cos(w*t)','y(0)=y0','Dy(0)=v0','t');

res=subs(y,{'y0','v0','w','theta','P','M','eps'},{1,2,w,theta,P,M,eps});

ezplot(char(res),[0 80 -1.5 1.5 ])

title ('analytic solution')

% подфункция вычисления правых частей уравнений

function F = oscil(t, y,theta,eps,b,w,P,M)

F = [y(2); -b*theta^2*sign(y(2))-2*eps*y(2)-theta^2*y(1)+P/M*cos(w*t)];

В результате получаем 2 графика: численного решения

и приближенного аналитического

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

Задание 4'

Задание 4.1' Построить графики свободных колебаний при заданных в таблице 1 значениях параметров при различных сочетаниях сухого и вязкого трения.

В качестве теста используем

, y(0)=1,

его решение y=exp(-t) (cos(2t)+3 sin(2t)).строим график

>> ezplot(@(t) exp(-t).*(cos(2*t)+3*sin(2*t)),[0 10 -0.5 2])

freeosc.mdl а вот и результат

как потом выяснилось)) можно задать и так..

4.2. Построить графики вынужденных колебаний системы с параметрами, заданными в таблице1.Рассмотреть случай гармонического воздействия и воздействия импульсного типа.

воздействия импульсного типа надо задавать как в методичке- с помощью SignalBuilder…лень ваще))

реализуем только гармоническое воздействие..пусть амплитуда возм. силы- 10…но в уравнении- 2*у….поэтому в окне Scope амплитуда установившихся колебаний- 2.

как-то так…forced_oscilation.mdl

4.3. Оформить скрипт-файл, содержащий результаты, полученные при выполнении заданий 1-3 и задания 4.Сопоставить и сравнить результаты.

это уж без меня))

Дата добавления: 2015-07-19