Шаматульский Роман P3212 Лаб 1-6

Р3212/shamatulskiy_409870/lab1/lab1.py

@@ -0,0 +1,69 @@
+import numpy as np
+
+def gauss_seidel(A, b, tol, max_iter=10000):
+    n = len(b)
+    x = np.zeros(n)
+    for k in range(1, max_iter+1):
+        x_old = x.copy()
+        for i in range(n):
+            s = np.dot(A[i, :i], x[:i]) + np.dot(A[i, i+1:], x_old[i+1:])
+            x[i] = (b[i] - s)/A[i, i]
+        err = np.abs(x - x_old)
+        if np.all(err < tol):
+            return x, k, err
+    return x, k, err
+
+def make_diagonally_dominant(A, b):
+    n = A.shape[0]
+    for i in range(n):
+        for j in range(i, n):
+            if abs(A[j, i]) >= np.sum(np.abs(A[j, :])) - abs(A[j, i]):
+                if j != i:
+                    A[[i, j]], b[[i, j]] = A[[j, i]].copy(), b[[j, i]].copy()
+                break
+    for i in range(n):
+        if abs(A[i, i]) < np.sum(np.abs(A[i, :])) - abs(A[i, i]):
+            return False
+    return True
+
+def read_matrix():
+    mode = input("Чтение матрицы из файла? (y/n): ")
+    if mode.lower() == "y":
+        fn = input("Имя файла: ")
+        data = np.loadtxt(fn)
+        A = data[:, :-1]
+        b = data[:, -1]
+    else:
+        n = int(input("n: "))
+        A = np.zeros((n, n))
+        b = np.zeros(n)
+        for i in range(n):
+            row = list(map(float, input(f"Строка {i+1} (n+1 чисел): ").split()))
+            A[i, :] = row[:-1]
+            b[i] = row[-1]
+    return A, b
+
+def read_tol():
+    mode = input("Чтение точности из файла? (y/n): ")
+    if mode.lower() == "y":
+        fn = input("Имя файла с точностью: ")
+        tol = float(open(fn).read().strip())
+    else:
+        tol = float(input("Точность: "))
+    return tol
+
+def main():
+    A, b = read_matrix()
+    tol = read_tol()
+    if not make_diagonally_dominant(A, b):
+        print("Невозможно добиться диагонального преобладания")
+        return
+    x, iters, err = gauss_seidel(A, b, tol)
+    norm = np.linalg.norm(A, ord=np.inf)
+    print("Норма матрицы (сумма по строкам, бесконечная):", norm)
+    print("Решение x:", x)
+    print("Число итераций:", iters)
+    print("Вектор погрешностей:", err)
+
+if __name__ == "__main__":
+    main()

Р3212/shamatulskiy_409870/lab2/lab2.py

@@ -0,0 +1,125 @@
+import math
+import sys
+import matplotlib.pyplot as plt
+
+def f1(x):
+    return x**3 - 2*x - 5
+
+def f2(x):
+    return math.sin(x) - 0.5*x
+
+def f3(x):
+    return math.exp(-x) - x
+
+functions = {
+    '1': ('x^3 - 2x - 5', f1),
+    '2': ('sin(x) - 0.5x', f2),
+    '3': ('e^(-x) - x', f3),
+}
+
+def phi3(x):
+    return math.exp(-x)
+
+phis = {
+    '3': phi3,
+}
+
+def bisection(f,a,b,eps):
+    fa,fb=f(a),f(b)
+    if fa*fb>0: return None,None,None
+    n=0
+    while (b-a)/2>eps:
+        c=(a+b)/2
+        if f(a)*f(c)<=0:
+            b=c
+        else:
+            a=c
+        n+=1
+    return (a+b)/2,f((a+b)/2),n
+
+def secant(f,a,b,eps):
+    fa,fb=f(a),f(b)
+    if fa==fb: return None,None,None
+    n=0
+    while abs(b-a)>eps:
+        c=b - fb*(b-a)/(fb-fa)
+        a,b,fa,fb=b,c,fb,f(c)
+        n+=1
+    return b,f(b),n
+
+def simple_iter(f,phi,a,b,eps):
+    mp = max(abs((phi(a+1e-6)-phi(a))/1e-6),abs((phi(b)-phi(b-1e-6))/1e-6))
+    if mp>=1: return 'bad',None,None
+    x0=(a+b)/2
+    n=0
+    while True:
+        x1=phi(x0)
+        n+=1
+        if abs(x1-x0)<eps:
+            return x1,f(x1),n
+        x0=x1
+
+def get_input():
+    if input('Ввод из файла? (y/n) ').lower()=='y':
+        fn=input('Имя файла: ')
+        with open(fn) as g:
+            data=g.read().split()
+        it=iter(data)
+        eq=next(it)
+        method=next(it)
+        a,b,eps=map(float,(next(it),next(it),next(it)))
+    else:
+        print('Уравнения:')
+        for k,v in functions.items():
+            print(k,v[0])
+        eq=input('Выберите номер функции: ')
+        print('Методы: 1-бисекция, 2-секущие, 3-простая итерация')
+        method=input('Выберите метод: ')
+        a=float(input('a= '))
+        b=float(input('b= '))
+        eps=float(input('ε= '))
+    return eq,method,a,b,eps
+
+def output(root,val,iter_count):
+    if input('Вывод в файл? (y/n) ').lower()=='y':
+        fn=input('Имя файла: ')
+        with open(fn,'w') as g:
+            g.write(f'Корень: {root}\n')
+            g.write(f'f(root): {val}\n')
+            g.write(f'Итераций: {iter_count}\n')
+    else:
+        print('Корень:',root)
+        print('f(root)=',val)
+        print('Итераций:',iter_count)
+
+def plot(f,a,b):
+    xs=[a+i*(b-a)/1000 for i in range(1001)]
+    ys=[f(x) for x in xs]
+    plt.plot(xs,ys)
+    plt.axhline(0, color='black')
+    plt.show()
+
+def main():
+    eq,method,a,b,eps=get_input()
+    desc,f=functions.get(eq,(None,None))
+    if f is None:
+        print('Неизвестная функция'); sys.exit()
+    fa,fb=f(a),f(b)
+    if method in ('1','2') and fa*fb>0:
+        print('На интервале нет единственного корня'); sys.exit()
+    if method=='1':
+        root,val,it=bisection(f,a,b,eps)
+    elif method=='2':
+        root,val,it=secant(f,a,b,eps)
+    else:
+        phi=phis.get(eq)
+        if not phi:
+            print('Итерационная форма не задана'); sys.exit()
+        root,val,it=simple_iter(f,phi,a,b,eps)
+        if root=='bad':
+            print('Условие сходимости не выполнено'); sys.exit()
+    output(root,val,it)
+    plot(f,a,b)
+
+if __name__=='__main__':
+    main()

Р3212/shamatulskiy_409870/lab2/lab2_system.py

@@ -0,0 +1,93 @@
+import math
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+def sys1(x):
+    u,v = x
+    return np.array([u**2 + v**2 - 4, u*v - 1])
+
+def phi1(x):
+    u,v = x
+    return np.array([math.sqrt(max(0,4 - v**2)), 1/u])
+
+def sys2(x):
+    u,v = x
+    return np.array([math.sin(u) + v - 1, u + math.cos(v) - 1])
+
+def phi2(x):
+    u,v = x
+    return np.array([math.asin(1 - v), math.acos(1 - u)])
+
+systems = {
+    '1': (sys1, phi1),
+    '2': (sys2, phi2),
+}
+
+def jacobian(phi, x, h=1e-6):
+    n = len(x)
+    J = np.zeros((n,n))
+    for i in range(n):
+        dx = np.zeros(n); dx[i] = h
+        J[:,i] = (phi(x+dx) - phi(x-dx))/(2*h)
+    return J
+
+def simple_iter(sysf, phi, x0, eps, maxit=1000):
+    J = jacobian(phi, x0)
+    if max(abs(np.linalg.eigvals(J))) >= 1:
+        print('Условие сходимости не выполнено'); sys.exit()
+    x_prev = np.array(x0)
+    errors = []
+    for k in range(1, maxit+1):
+        x_new = phi(x_prev)
+        err = np.abs(x_new - x_prev)
+        errors.append(err)
+        if np.all(err < eps):
+            return x_new, k, errors
+        x_prev = x_new
+    print('Не сошлось за max итераций'); sys.exit()
+
+def get_input():
+    print('Системы:')
+    print('1: x^2 + y^2 = 4; x*y = 1')
+    print('2: sin(x) + y = 1; x + cos(y) = 1')
+    sysn = input('Выберите номер системы: ')
+    x0 = float(input('x0= ')), float(input('y0= '))
+    eps = float(input('ε= '))
+    return sysn, x0, eps
+
+def output(x, it, err):
+    if input('Вывод в файл? (y/n) ').lower()=='y':
+        fn = input('Имя файла: ')
+        with open(fn, 'w') as g:
+            g.write(f'x = {x}\n')
+            g.write(f'Итераций = {it}\n')
+            g.write(f'Ошибка последней итерации = {err}\n')
+    else:
+        print('x =', x)
+        print('Итераций =', it)
+        print('Ошибка последней итерации =', err)
+
+def plot(sysf, a=-2, b=2):
+    xs = np.linspace(a, b, 400)
+    ys = np.linspace(a, b, 400)
+    X, Y = np.meshgrid(xs, ys)
+    Z1 = sysf([X, Y])[0]
+    Z2 = sysf([X, Y])[1]
+    plt.contour(X, Y, Z1, levels=[0])
+    plt.contour(X, Y, Z2, levels=[0])
+    plt.axhline(0, color='black')
+    plt.axvline(0, color='black')
+    plt.show()
+
+def main():
+    sysn, x0, eps = get_input()
+    sysf, phi = systems.get(sysn, (None,None))
+    if not sysf:
+        print('Неизвестная система'); sys.exit()
+    x, it, errors = simple_iter(sysf, phi, x0, eps)
+    output(x, it, list(errors[-1]))
+    plot(sysf)
+
+if __name__=='__main__':
+    main()