info2.Rmd

---
title: "Układ ze sprężyną o nieliniowej charakterystyce"
author:
- name: Jan Adamski
  nr: 123456
- name: Adam Jański
  nr: 654321
output:
  pdf_document:
    keep_tex: yes
    template: spraw_templ.tex
  html_document:
    self_contained: no
  word_document: default
class: Informatyka II
deadline: 01.06.2019 r.
instructor: Jakiś Ktoś
---

```{r global_options, include=FALSE}
# Opcje uruchamiania R
knitr::opts_chunk$set(
  warning=FALSE, # Nie wyświetlaj ostrzeżeń
  echo=FALSE, # Powoduje nie umieszczanie kodu R w raporcie
  fig.align='center', # Wyśrodkowane wykresy
  fig.height = 4, # Wysokość wykresów
  fig.width = 5, # Szerokość wykresów
  dev.args=list(pointsize=10) # Rozmiar czcionki na wykresach
)
```

# Opis problemu
```{r scheme}
# Schemat masy na sprężynce
par(mar=c(4.1,0.5,0.5,0.5))
plot(NA,xlim=c(-3,3),ylim=c(-1,1),asp=1,xlab="X",ylab="",bty="n",yaxt="n")
rect(0,-0.5,1,0.5)
segments(-3,-1,-3,1)
rect(-3.3,-1,-3,1,density=10,border = NA)
segments(0.5,-0.5,0.5,-2,lty=2)
n=10
lines(cbind(seq(-3,0,len=2*n+4),c(0,0,rep(c(1,-1)/10,n),0,0)+0.2))
text(-1.5,0.5,"k")
lines(cbind(c(-3,-1.5,-1.5,-1.5),c(0,0,-0.15,0.15)-0.2))
lines(cbind(c(-1.6,-1.4,-1.4,0,-1.4,-1.4,-1.6),c(-0.2,-0.2,0,0,0,0.2,0.2)-0.2))
text(-1.5,-0.6,"c")
arrows(0.5,0.2,0.2,0.2,length = 0.1)
text(1.1,0.2,"Fs=-kx",adj=c(0,0.5))
arrows(0.5,-0.2,0.3,-0.2,length = 0.1)
text(1.1,-0.2,"Ft=-cv",adj=c(0,0.5))
arrows(0.5,0,0.8,0,length = 0.1)
text(1.1,0,"v",adj=c(0,0.5))
```

Charakterystyka nieliniowa sprężyny: $k=k_1(1+k_2x^2)$

# Równania ruchu

Siła wywierana na masę przez sprężynę:
\[F_s = -kx = -k_1(1+k_2x^2)x\]
Siła wywierana na masę przez tłumik:
\[F_t = -c\dot x\]
Równanie ruchu:
\[m\ddot x = F_s + F_t = -k_1(1+k_2x^2)x-c\dot x\]
Równanie to można przekształcić na układ równań pierwszego rzędu:
\begin{equation}\label{uklad}\begin{cases}
\dot x &= y\\
\dot y &= \frac{1}{m}\left(-k_1(1+k_2x^2)x - cy\right)
\end{cases}\end{equation}

Energia kinetyczna:
\[E_k = \frac{m\dot x^2}{2}\]
Energia potencjalna sprężyny:
\[E_p = \int k_1(1+k_2x^2)x\cdot dx = \int \left(k_1x+k_1k_2x^3\right)dx = \frac{k_1x^2}{2}+\frac{k_1k_2x^4}{4}\]
Całkowita energia mechaniczna:
\[E = \frac{m\dot x^2}{2} + \frac{k_1x^2}{2}+\frac{k_1k_2x^4}{4}\]

**Uwaga:** bez straty ogólności, w dalszej części pracy będziemy przyjmować $m=1$.

# Metoda Obliczeniowa

Układ równań został scałkowany przy pomocy metody Runge-Kutta 4-tego rzędu. Czas całkowania: $10s$. Krok całkowania: $\frac{10s}{1500}$.

# Wyniki

Przeprowadzono symulację dla czterech przypadków: Liniowej i nieliniowej charakterystyki sprężyny, a także z i bez tłumienia.

```{r}
# Wszystkie wykresy bez górnego i prawego marginesu
par(mar=c(4.1,4.1,0.5,0.5))
```

## Przypadek liniowej sprężyny ($k_1=10$)
```{r}
# Wczytanie danych wygenerowanych przez program
  k1=10
  k2=0
  c=0
  tab = read.csv("info2_program/linear.csv")
  tab$E = (tab$V^2)/2 + (k1*tab$X^2)/2+(k1*k2*tab$X^4)/4
```

Liniowa sprężyna charakteryzuje się liniową zależnością pomiędzy wychyleniem a siłą:

```{r linear-characteristic}
x=seq(-1.5,1.5,len=200)
plot(x,x*k1,type="l",xlab="x [m]",ylab="F [N]")
abline(h=0,v=0,lty=2)
```
Rozwiązanie numeryczne tego układu jest zgodne z oczekiwanym sinusoidalnym kształtem:

```{r linear-solution}
  plot(tab$T,tab$X,xlab="t [s]",ylab="x [m]",type="l",lwd=2)
  abline(h=0,lty=2)
```

Dla tak prostego układu $m\ddot x + k_1x = 0$, możemy wyznaczyć rozwiązanie analityczne przez podstawienie $x=e^{rt}$. Otrzymane rozwiązanie ogólne (po usunięciu części urojonej) ma postać:
\[x(t) = A \cos{\sqrt{\frac{k_1}{m}}t} + B \sin{\sqrt{\frac{k_1}{m}}t}\]
Dla warunków początkowych $x(0)=1\land \dot x(0)=1$ mamy:
\[x(t) = \cos{\sqrt{\frac{k_1}{m}}t}\]

Wizualnie rozwiązanie numeryczne pokryłoby się z analitycznym, możemy jednak zwizualizować błąd względny ($\left|\frac{x_n-x_a}{x_a}\right|$) w skali logarytmicznej:

```{r linear-error}
  tab$Xa=cos(sqrt(k1)*tab$T)
  plot(tab$T,abs((tab$X-tab$Xa)/tab$Xa),xlab="t [s]",ylab="Błąd względny [%]",type="l",lwd=2,log="y")
```

W przestrzeni fazowej ($x-v$) rozwiązanie jest zamkniętą elipsą. Nie ma tu sensu mówienie o jej proporcjach, ponieważ obie osie mają różne skale.

```{r linear-phase-space}
plot(tab$X, tab$V, xlab="x [m]", ylab="v [m/s]", type="l", lwd=2, xlim=c(-1,1), ylim=c(-4,4))
```

Zamknięta ścieżka w przestrzeni fazowej sugeruje, że układ nie traci energii.

## Przypadek nieliniowej sprężyny ($k_1=5\quad k_2=3$)
```{r}
  k1=5
  k2=3
  c=0
  tab = read.csv("info2_program/nonlinear.csv")
  tab$E = (tab$V^2)/2 + (k1*tab$X^2)/2+(k1*k2*tab$X^4)/4
```

Zależność siły od wychylenia dla omawianej nieliniowej sprężyny:

```{r nonlinear-characteristic}
x=seq(-1.5,1.5,len=200)
plot(x,x*(k1+k1*k2*x^2),type="l",xlab="x [m]",ylab="F [N]")
abline(h=0,v=0,lty=2)
abline(0,k1,lty=3)
```

Rozwiązanie numeryczne dla nieliniowej sprężyny ma "ostrzej" zakończone maxima i minima. Jest to związane z wyższą siłą siłą przy dużych wychyleniach niż w przypadku liniowej sprężyny.

```{r nonlinear-solution}
  plot(tab$T,tab$X,xlab="t [s]",ylab="x [m]",type="l",lwd=2)
  abline(h=0,lty=2)
```

W przestrzeni fazowej, trajektoria nadal jest zamkniętą pętlą, lecz nie jest już elipsą:

```{r nonlinear-phase}
plot(tab$X, tab$V, xlab="x [m]", ylab="v [m/s]", type="l", lwd=2, xlim=c(-1,1), ylim=c(-4,4))
```

## Przypadek liniowej sprężyny z tłumikiem ($k_1=10\quad c=0.3$)

```{r}
  k1=10
  k2=0
  c=0.3
  tab = read.csv("info2_program/linear_dumper.csv")
  tab$E = (tab$V^2)/2 + (k1*tab$X^2)/2+(k1*k2*tab$X^4)/4
```

Rozwiązanie numeryczne dla liniowej sprężyny z tłumikiem, ma wykładniczy spadek. W związku z tym, że częstotliwość drgania układu liniowego nie zależy od wychylenia, odstępy pomiędzy momentami przejścia przez zero są ciągle stałe.

```{r linear-dump-solution}
  plot(tab$T,tab$X,xlab="t [s]",ylab="x [m]",type="l",lwd=2)
  abline(h=0,lty=2)
```

W przestrzeni fazowej, trajektoria nie jest już zamknięta i schodzi spiralnie do zera:

```{r linear-dump-phase}
plot(tab$X, tab$V, xlab="x [m]", ylab="v [m/s]", type="l", lwd=2, xlim=c(-1,1), ylim=c(-4,4))
```

Tłumik wprowadza element niezachowawczy do układu. Powoduje to spadek energii. Łatwo ten efekt zobaczyć wymnażając oryginalny układ przez $\dot x$:
\[\dot xm\ddot x + \dot xc\dot x + \dot xk_1x=0\]
Po przekształceniu mamy:
\[\frac{d}{dt}\left(\frac{m\dot x^2}{2}+\frac{k_1 x^2}{2}\right)=-c\dot x^2\]
Na wykresie całkowitej energii widać wyraźnie spadek energii, który następuje w momentach wysokiej prędkości.

```{r linear-dump-energy}
plot(tab$T, tab$E, xlab="x [m]", ylab="E [J]", type="l", lwd=2, ylim=c(0,7))
```

## Przypadek nieliniowej sprężyny z tłumikiem ($k_1=5\quad k_2=3\quad c=0.3$)

```{r}
  k1=5
  k2=3
  c=0.3
  tab = read.csv("info2_program/nonlinear_dumper.csv")
  tab$E = (tab$V^2)/2 + (k1*tab$X^2)/2+(k1*k2*tab$X^4)/4
```

Rozwiązanie numeryczne dla nie-liniowej sprężyny z tłumikiem, ma także wykładniczy spadek. Dla nieliniowej sprężyny częstotliwość zmienia się wraz z maksymalnym wychyleniem. Dlatego odstępy pomiędzy przejściami przez zero będą się wydłużać.

```{r nonlinear-dump-solution}
  plot(tab$T,tab$X,xlab="t [s]",ylab="x [m]",type="l",lwd=2)
  abline(h=0,lty=2)
```

W przestrzeni fazowej, trajektoria jest nieeliptyczną spiralą. Dodatkowo można zauważyć że gdy wychylenie się zmniejsza, trajektoria robi się coraz bardziej eliptyczna:

```{r nonlinear-dump-phase}
plot(tab$X, tab$V, xlab="x [m]", ylab="v [m/s]", type="l", lwd=2, xlim=c(-1,1), ylim=c(-4,4))
```

Tak jak w poprzednim przypadku, energia spada w czasie:

```{r nonlinear-dump-energy}
plot(tab$T, tab$E, xlab="x [m]", ylab="E [J]", type="l", lwd=2, ylim=c(0,7))
```

# Omówienie wyników