From a163280f7a0a82335b0705ecfb935a18f980bfd3 Mon Sep 17 00:00:00 2001 From: Rokas Puzonas Date: Thu, 21 Dec 2023 15:50:48 +0200 Subject: [PATCH] solve lab4 --- Lab3/part_1.py | 22 ++--- Lab3/part_2.py | 6 +- Lab3/part_4.py | 2 +- Lab4/main.py | 227 +++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 239 insertions(+), 18 deletions(-) create mode 100644 Lab4/main.py diff --git a/Lab3/part_1.py b/Lab3/part_1.py index 74244b9..c898a54 100644 --- a/Lab3/part_1.py +++ b/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) - #plot_vienanariu_paklaida(F, x_s, vienanariai) + #plt.plot(mazgai[0], mazgai[1], "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))) - 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 ) \ No newline at end of file diff --git a/Lab3/part_2.py b/Lab3/part_2.py index b9ff205..260acf2 100644 --- a/Lab3/part_2.py +++ b/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): diff --git a/Lab3/part_4.py b/Lab3/part_4.py index 37b8b20..b8c920f 100644 --- a/Lab3/part_4.py +++ b/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, ) \ No newline at end of file diff --git a/Lab4/main.py b/Lab4/main.py new file mode 100644 index 0000000..1d4527d --- /dev/null +++ b/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 +) \ No newline at end of file