
clear;

# ======================================================================
# globale Variablen
# ======================================================================

global k
global phi_0 
global t_0 

# ======================================================================
# Basisfunktion Psi_0(t)
# ======================================================================

function value = psi0_0 (x)

  global phi_0;
  global t_0;

  value = phi_0 / t_0 * x;
    
endfunction

##function value = psi0_n (n, x)
##
##  global a;
##  value = x/a - (x/a).^(n+1);
##    
##endfunction
##
##function value = psi2_n (n, x)
##  
##  global a;
##  value = - (n * (n+1)) / (a^(n+1)) * x.^(n-1);
##    
##endfunction

# ======================================================================
# Basisfunktion Psi_n(t)
# ======================================================================

function value = psi0_n (n, x)

  global t_0;
  value = (x / t_0) - (x / t_0).^(2 * n + 1);
    
endfunction

# ======================================================================
# Zweite Abgleiutung der Basisfunktion Psi_n(t)
# ======================================================================

function value = psi2_n (n, x)
  
  global t_0;
  value = - (2 * n + 1) * (2 * n) / (t_0^(2 * n + 1)) * x.^(2 * n - 1);
    
endfunction

# ======================================================================
# h(x) = g(t) - L{Psi_0(t)}  mit g(t) = 0 für diesen Fall
# ======================================================================

function value = h(x)

  global k;
  global t_0;
  global phi_0;
  
  value = -k * phi_0 / t_0 * x;
  
endfunction

# ======================================================================
# g_n(x) = L{psi_n(x)}
# ======================================================================

function value = g_n(n, x)

  global k;
  value = k * psi0_n(n, x) + psi2_n(n, x);
  
endfunction


# ======================================================================
# f(x) = psi_0(x) + a_1 * psi_1(x) + a_2 * psi_2(x) + ...
# ======================================================================

function value = f(alpha, x)
  
  value = psi0_0(x);
  
  for n = 1:length(alpha)
    value = value + alpha(n) * psi0_n(n, x);
  endfor  
    
endfunction

# ======================================================================
# Matrix aufstellen
# ======================================================================

function M = Matrix(x)

  N = length(x)

  for n = 1:N
    for j = 1:N
     
      M(j, n) = g_n(n, x(j));
    
    endfor
  endfor

endfunction

# ======================================================================
# Plot erstellen
# ======================================================================

function makePlot(x1, y1, x2, y2)
  
  hold off;
  figure(1, 'position',[100,100,500,300]);
  plot(x1,y1);
  hold on;
  grid on;
  plot(x2,y2);
  
  xlabel('t/s');
  ylabel('\phi(t)');
 
endfunction

# ======================================================================
# Lösung der Differentialgleichung für das Fadenpendel
# ======================================================================

# Vorgaben
# ----------------------------------------------------------------------

g     = 9.81;                     # Erdbeschleunigung
l     = 2;                        # Länge des Fadenpendels

T = 2 * pi * sqrt(l / g);         # Schwingungsdauer

t_0     = T / 4;                  # beliebiger Zeitpunkt (hier T/4)
phi_0   = 1;                      # Auslenkung des Pendels bei t_0


k = g/l                           # Hilfskonstante um Schreibarbeit zu sparen

# Lösung der Differentialgleichung mit der Momentenmethode
# ----------------------------------------------------------------------

x   = linspace(0.1 * t_0, t_0, 3)'; # Testpunkte erstellen

M   = Matrix(x)                     # Matrix für Testpunkte aufstellen
vh  = h(x)

va  = inv(M) * vh                   # Gleichungssystem lösen (alpha-Werte)

px  = linspace(-2 * t_0, 2 * t_0, 50);

py_1 = f(va, px);                   # Nährungslösung durch Momentenmethode
py_2 = phi_0 * sin (2 * pi / T * px);  # exakte Lösung

makePlot(px, py_1, px, py_2);

# ----------------------------------------------------------------------



