solve lab4

Lab3/part_1.py

@@ -44,10 +44,6 @@ def gauso_atispindzio_metodas(lygciu_sistema: LygciuSistema, epsilon=1e-7):
         A1[i:n, :] = Q.dot(A1[i:n, :])
 
     if np.sum(np.abs(A1[n-1, 0:n-nb+1])) < epsilon:
-        # if abs(A1[n-1, n]) < epsilon:
-        #     print("Be galo daug sprendiniu")
-        # else:
-        #     print("Nera sprendiniu")
         return None
 
     # atgalinis etapas:
@@ -55,12 +51,6 @@ def gauso_atispindzio_metodas(lygciu_sistema: LygciuSistema, epsilon=1e-7):
     for i in range(n - 1, -1, -1):
         x[i, :] = (A1[i, n:n + nb] - A1[i, i + 1:n] * x[i + 1:n, :]) / A1[i, i]
 
-    # print("Sprendinys:")
-    # for i in range(0, n):
-    #     print(f"x{i} = {x[i, 0]:.5f}")
-
-    # print("Sklaida:", gauti_sklaida(lygciu_sistema, list(x.flat)))
-
     return list(x.flat)
 
 def gauti_x_taskus(from_x, to_x, density) -> list[float]:
@@ -142,22 +132,24 @@ def plot_vienanariu_paklaida(F, x_s: list[float], vienanariai):
 def main(F, x_range, mazgu_kiekis, tarpiniai_taskai):
     from_x, to_x = x_range
 
-    if True:
+    if False:
         mazgai = gauti_xy_taskus(F, from_x, to_x, mazgu_kiekis)
         vienanariai = gauti_vienanario_sprendinius(mazgai[0], mazgai[1])
         x_s = gauti_x_taskus(from_x, to_x, tarpiniai_taskai)
 
-        plot_vienanariu_palyginima(F, x_s, vienanariai)
+        #plt.plot(mazgai[0], mazgai[1], "o", label="Mazgas")
-        #plot_vienanariu_paklaida(F, x_s, vienanariai)
+        #plot_vienanariu_palyginima(F, x_s, vienanariai)
+        plot_vienanariu_paklaida(F, x_s, vienanariai)
         print("Paklaida:", sum(gauti_vienanariu_paklaidas(F, x_s, vienanariai)))
 
-    if False:
+    if True:
         ciobysevo_abscises = gauti_ciobysevo_abscises(mazgu_kiekis)
         x_mazgai = konvertuoti_is_ciobysevo(ciobysevo_abscises, to_x, from_x)
         y_mazgai = list(F(x) for x in x_mazgai)
         vienanariai = gauti_vienanario_sprendinius(x_mazgai, y_mazgai)
 
         x_s = gauti_x_taskus(from_x, to_x, tarpiniai_taskai)
+        plt.plot(x_mazgai, y_mazgai, "o", label="Mazgas")
         plot_vienanariu_palyginima(F, x_s, vienanariai)
         #plot_vienanariu_paklaida(F, x_s, vienanariai)
         print("Paklaida:", sum(gauti_vienanariu_paklaidas(F, x_s, vienanariai)))
@@ -168,5 +160,5 @@ main(
     F = lambda x: cos(2*x) * (sin(3*x) + 1.5) - cos(x/5),
     x_range = (-2, 3),
     mazgu_kiekis = 10,
-    tarpiniai_taskai = 30
+    tarpiniai_taskai = 200
 )

Lab3/part_2.py

@@ -74,7 +74,8 @@ def draw_hermite_curve(X: list[float], Y: list[float], dY: list[float], scalar =
     assert len(X) == len(Y) == len(dY)
     N = len(X)
 
-    plt.plot(X[0], Y[0], 'ro')
+    plt.plot(X, Y, '--g', label="Originali")
+    plt.plot(X[0], Y[0], 'ro', label="Mazgas")
     for i in range(N - 1):
         plot_x = np.linspace(X[i], X[i+1], scalar)
         plot_y = []
@@ -91,8 +92,9 @@ def draw_hermite_curve(X: list[float], Y: list[float], dY: list[float], scalar =
             f  = U1*Y[i  ] + V1*dY[i  ]
             f += U2*Y[i+1] + V2*dY[i+1]
             plot_y.append(f)
-        plt.plot(plot_x, plot_y, 'b')
+        plt.plot(plot_x, plot_y, 'b', label="Interpoliuota" if i == 0 else None)
         plt.draw()
+    plt.legend()
     plt.show()
 
 def main(country, interpolation):

Lab3/part_4.py

@@ -139,5 +139,5 @@ def main(country: str, levels: list[int], show_accuracy_plot):
 main(
     country = "Zambia",
     levels = [1, 4, 8, 10],
-    show_accuracy_plot = False,
+    show_accuracy_plot = True,
 )

Lab4/main.py

@@ -0,0 +1,227 @@
+from typing import Literal
+import numpy as np
+import matplotlib.pyplot as plt
+from dataclasses import dataclass
+import scipy.integrate
+
+@dataclass
+class Salyga:
+    m1: float # Parašiutininko masė
+    m2: float # Įrangos masė
+    h0: float # Iššokimo aukštis
+    tg: float # Laikas iki parašiuto išskleidimo
+    k1: float # Oro pasipriešinimas be parašiuto
+    k2: float # Oro pasipriešinimas su parašiutu
+    g = 9.81
+
+def ispresti_euleriu(salyga: Salyga, iteracijos: float, simuliacijos_laikotarpis: float):
+    t_history = []
+    h_history = []
+    v_history = []
+
+    m = salyga.m1 + salyga.m2
+    dt = simuliacijos_laikotarpis / iteracijos
+    h = salyga.h0
+    t = 0
+    v = 0
+    for _ in range(iteracijos):
+        k = salyga.k1 if t < salyga.tg else salyga.k2
+        pagreitis = salyga.g - k * v ** 2 / m
+
+        h += v * dt
+        v += -pagreitis * dt
+
+        if h <= 0:
+            h = 0
+            v = 0
+
+        t_history.append(t)
+        h_history.append(h)
+        v_history.append(v)
+
+        t += dt
+
+    return t_history, h_history, v_history
+
+def ispresti_rk4(salyga: Salyga, iteracijos: int, simuliacijos_laikotarpis: float):
+    def funk(X, t):
+        nonlocal salyga
+        k = salyga.k1 if t < salyga.tg else salyga.k2
+        v = X[1]
+        m = salyga.m1 + salyga.m2
+        pagreitis = -salyga.g + k * v ** 2 / m
+
+        return np.array([v, pagreitis])
+
+    t = np.linspace(0, simuliacijos_laikotarpis, iteracijos)
+    dt = t[1]-t[0]
+    rez = np.zeros([2, iteracijos], dtype=float)
+    rez[:,0] = np.array([salyga.h0, 0])
+
+    for i in range(iteracijos-1):
+        fz   = rez[:,i] + funk(rez[:,i],t[i]     ) * dt/2
+        fzz  = rez[:,i] + funk(fz      ,t[i]+dt/2) * dt/2
+        fzzz = rez[:,i] + funk(fzz     ,t[i]+dt/2) * dt
+
+        rez[:,i+1] = rez[:,i] + dt/6 * (
+                funk(rez[:,i],t[i]     ) +
+            2 * funk(fz      ,t[i]+dt/2) +
+            2 * funk(fzz     ,t[i]+dt/2) +
+                funk(fzzz    ,t[i]+dt  )
+        )
+
+        if rez[0,i+1] <= 0:
+            rez[0,i+1] = 0
+            rez[1,i+1] = 0
+
+    return t, rez[0,:], rez[1,:]
+
+def main_1(salyga: Salyga, iteracijos: int, simuliacijos_laikotarpis):
+    t_history_e  , h_history_e  , v_history_e   = ispresti_euleriu(salyga, iteracijos, simuliacijos_laikotarpis)
+    t_history_rk4, h_history_rk4, v_history_rk4 = ispresti_rk4(salyga, iteracijos, simuliacijos_laikotarpis)
+
+    print("Žingsnio dydis: ", simuliacijos_laikotarpis / iteracijos)
+
+    for i, h in enumerate(h_history_e):
+        if h == 0:
+            print("Eulerio", v_history_e[i-1], t_history_e[i-1])
+            break
+
+    for i, h in enumerate(h_history_rk4):
+        if h == 0:
+            print("rk4", v_history_rk4[i-1], t_history_rk4[i-1])
+            break
+
+    fig1=plt.figure(1)
+
+    ax1 = fig1.add_subplot(1,2,1)
+    #ax1.plot(t_history_e, h_history_e, 'r-', label="Eulerio")
+    ax1.plot(t_history_rk4, h_history_rk4, 'b-', label="IV eilės Rungės ir Kutos")
+    #ax1.legend()
+    ax1.set_xlabel("t (s)")
+    ax1.set_ylabel("h (m)")
+    ax1.set_title("Aukštis")
+
+    ax2 = fig1.add_subplot(1,2,2)
+    #ax2.plot(t_history_e, v_history_e, 'r-', label="Eulerio")
+    ax2.plot(t_history_rk4, v_history_rk4, 'b-', label="IV eilės Rungės ir Kutos")
+    #ax2.legend()
+    ax2.set_xlabel("t (s)")
+    ax2.set_ylabel("v (m/s)")
+    ax2.set_title("Greitis")
+
+    plt.show()
+
+def main_2(salyga: Salyga, metodas: Literal["euler", "rk4"], iteracijos: list[int], simuliacijos_laikotarpis):
+    fig1=plt.figure(1)
+
+    ax1 = fig1.add_subplot(1,2,1)
+    ax1.set_xlabel("t (s)")
+    ax1.set_ylabel("h (m)")
+    ax1.set_title("Aukštis")
+
+    ax2 = fig1.add_subplot(1,2,2)
+    ax2.set_xlabel("t (s)")
+    ax2.set_ylabel("v (m/s)")
+    ax2.set_title("Greitis")
+
+    cmap = plt.cm.get_cmap('hsv', len(iteracijos)+1)
+    for i, iteraciju_kiekis in enumerate(iteracijos):
+        if metodas == "euler":
+            t_history, h_history, v_history = ispresti_euleriu(salyga, iteraciju_kiekis, simuliacijos_laikotarpis)
+        elif metodas == "rk4":
+            t_history, h_history, v_history = ispresti_rk4(salyga, iteraciju_kiekis, simuliacijos_laikotarpis)
+        zingsnis = simuliacijos_laikotarpis / iteraciju_kiekis
+        ax1.plot(t_history, h_history, c=cmap(i))
+        ax2.plot(t_history, v_history, c=cmap(i), label=f"{zingsnis:.4f}")
+
+    fig1.legend()
+    plt.show()
+
+def main_3(salyga: Salyga, iteracijos, simuliacijos_laikotarpis: float):
+    print("Žingsnio dydis: ", simuliacijos_laikotarpis / iteracijos)
+
+    t_history_e  , h_history_e  , v_history_e   = ispresti_euleriu(salyga, iteracijos, simuliacijos_laikotarpis)
+    t_history_rk4, h_history_rk4, v_history_rk4 = ispresti_rk4(salyga, iteracijos, simuliacijos_laikotarpis)
+
+    def funk(t, X):
+        nonlocal salyga
+        k = salyga.k1 if t < salyga.tg else salyga.k2
+        v = X[1]
+        m = salyga.m1 + salyga.m2
+        pagreitis = -salyga.g + k * v ** 2 / m
+
+        return np.array([v, pagreitis])
+
+    tspan = np.array([0, simuliacijos_laikotarpis])
+    Y = scipy.integrate.solve_ivp(funk, tspan, [salyga.h0, 0])
+
+    fig1=plt.figure(1)
+
+    zero_point = 0
+    for i, h in enumerate(Y.y[0,:]):
+        if h <= 0:
+            zero_point = i
+            break
+
+    ax1 = fig1.add_subplot(1,2,1)
+    ax1.set_xlabel("t (s)")
+    ax1.set_ylabel("h (m)")
+    ax1.set_title("Aukštis")
+    ax1.plot(Y.t, Y.y[0,:], color="r", label="scipy")
+    ax1.plot(t_history_e, h_history_e, color="g", label="Eulerio")
+    ax1.plot(t_history_rk4, h_history_rk4, color="b", label="IV eilės Rungės ir Kutos")
+    ax1.plot(Y.t[zero_point], Y.y[0, zero_point], 'k.')
+    ax1.legend()
+
+    ax2 = fig1.add_subplot(1,2,2)
+    ax2.set_xlabel("t (s)")
+    ax2.set_ylabel("v (m/s)")
+    ax2.set_title("Greitis")
+    ax2.plot(Y.t, Y.y[1,:], color="r", label="scipy")
+    ax2.plot(t_history_e, v_history_e, color="g", label="Eulerio")
+    ax2.plot(t_history_rk4, v_history_rk4, color="b", label="IV eilės Rungės ir Kutos")
+    ax2.plot(Y.t[zero_point], Y.y[1,zero_point], 'k.')
+    ax2.legend()
+
+    plt.show()
+
+# Variantas 20
+salyga = Salyga(
+    m1 = 120.0,
+    m2 = 15.0,
+    h0 = 2800.0,
+    tg = 30.0,
+    k1 = 0.15,
+    k2 = 10.0
+)
+
+main_1(
+    salyga,
+    iteracijos = 20000,
+    simuliacijos_laikotarpis = 100
+)
+
+main_2(
+    salyga,
+
+    # Žemi žingsniai
+    # metodas="rk4",
+    # iteracijos = [1000, 10000, 20000, 40000],
+    # metodas="euler",
+    # iteracijos = [1000, 10000, 20000, 40000],
+
+    # Aukšti žingsniai
+    metodas="rk4",
+    iteracijos = [103, 120, 200, 500, 1000, 10000, 40000],
+    # metodas="euler",
+    # iteracijos = [609, 700, 1000, 10000, 20000],
+
+    simuliacijos_laikotarpis = 100
+)
+
+main_3(
+    salyga,
+    iteracijos = 20000,
+    simuliacijos_laikotarpis = 100
+)