diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b694934 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.venv \ No newline at end of file diff --git a/Pipfile b/Pipfile index b9ba84f..ecfc0ab 100644 --- a/Pipfile +++ b/Pipfile @@ -4,6 +4,11 @@ verify_ssl = true name = "pypi" [packages] +"pyqt5" = "*" +sympy = "*" +numpy = "*" +scipy = "*" +matplotlib = "*" [dev-packages] diff --git a/Pipfile.lock b/Pipfile.lock new file mode 100644 index 0000000..48cebb6 --- /dev/null +++ b/Pipfile.lock @@ -0,0 +1,205 @@ +{ + "_meta": { + "hash": { + "sha256": "819dbb4221941078ef4ee76f8d1baabd283dc82604457655f6cdb74f4561f22c" + }, + "pipfile-spec": 6, + "requires": { + "python_version": "3.7" + }, + "sources": [ + { + "name": "pypi", + "url": "https://pypi.org/simple", + "verify_ssl": true + } + ] + }, + "default": { + "cycler": { + "hashes": [ + "sha256:1d8a5ae1ff6c5cf9b93e8811e581232ad8920aeec647c37316ceac982b08cb2d", + "sha256:cd7b2d1018258d7247a71425e9f26463dfb444d411c39569972f4ce586b0c9d8" + ], + "version": "==0.10.0" + }, + "kiwisolver": { + "hashes": [ + "sha256:0ee4ed8b3ae8f5f712b0aa9ebd2858b5b232f1b9a96b0943dceb34df2a223bc3", + "sha256:0f7f532f3c94e99545a29f4c3f05637f4d2713e7fd91b4dd8abfc18340b86cd5", + "sha256:1a078f5dd7e99317098f0e0d490257fd0349d79363e8c923d5bb76428f318421", + "sha256:1aa0b55a0eb1bd3fa82e704f44fb8f16e26702af1a073cc5030eea399e617b56", + "sha256:2874060b91e131ceeff00574b7c2140749c9355817a4ed498e82a4ffa308ecbc", + "sha256:379d97783ba8d2934d52221c833407f20ca287b36d949b4bba6c75274bcf6363", + "sha256:3b791ddf2aefc56382aadc26ea5b352e86a2921e4e85c31c1f770f527eb06ce4", + "sha256:4329008a167fac233e398e8a600d1b91539dc33c5a3eadee84c0d4b04d4494fa", + "sha256:45813e0873bbb679334a161b28cb9606d9665e70561fd6caa8863e279b5e464b", + "sha256:53a5b27e6b5717bdc0125338a822605084054c80f382051fb945d2c0e6899a20", + "sha256:574f24b9805cb1c72d02b9f7749aa0cc0b81aa82571be5201aa1453190390ae5", + "sha256:66f82819ff47fa67a11540da96966fb9245504b7f496034f534b81cacf333861", + "sha256:79e5fe3ccd5144ae80777e12973027bd2f4f5e3ae8eb286cabe787bed9780138", + "sha256:83410258eb886f3456714eea4d4304db3a1fc8624623fc3f38a487ab36c0f653", + "sha256:8b6a7b596ce1d2a6d93c3562f1178ebd3b7bb445b3b0dd33b09f9255e312a965", + "sha256:9576cb63897fbfa69df60f994082c3f4b8e6adb49cccb60efb2a80a208e6f996", + "sha256:95a25d9f3449046ecbe9065be8f8380c03c56081bc5d41fe0fb964aaa30b2195", + "sha256:a424f048bebc4476620e77f3e4d1f282920cef9bc376ba16d0b8fe97eec87cde", + "sha256:aaec1cfd94f4f3e9a25e144d5b0ed1eb8a9596ec36d7318a504d813412563a85", + "sha256:acb673eecbae089ea3be3dcf75bfe45fc8d4dcdc951e27d8691887963cf421c7", + "sha256:b15bc8d2c2848a4a7c04f76c9b3dc3561e95d4dabc6b4f24bfabe5fd81a0b14f", + "sha256:b1c240d565e977d80c0083404c01e4d59c5772c977fae2c483f100567f50847b", + "sha256:c595693de998461bcd49b8d20568c8870b3209b8ea323b2a7b0ea86d85864694", + "sha256:ce3be5d520b4d2c3e5eeb4cd2ef62b9b9ab8ac6b6fedbaa0e39cdb6f50644278", + "sha256:e0f910f84b35c36a3513b96d816e6442ae138862257ae18a0019d2fc67b041dc", + "sha256:ea36e19ac0a483eea239320aef0bd40702404ff8c7e42179a2d9d36c5afcb55c", + "sha256:efabbcd4f406b532206b8801058c8bab9e79645b9880329253ae3322b7b02cd5", + "sha256:f923406e6b32c86309261b8195e24e18b6a8801df0cfc7814ac44017bfcb3939" + ], + "version": "==1.0.1" + }, + "matplotlib": { + "hashes": [ + "sha256:1ae6549976b6ceb6ee426272a28c0fc9715b3e3669694d560c8f661c5b39e2c5", + "sha256:4d4250bf508dd07cca3b43888097f873cadb66eec6ac63dbbfb798798ec07af2", + "sha256:53af2e01d7f1700ed2b64a9091bc865360c9c4032f625451c4589a826854c787", + "sha256:63e498067d32d627111cd1162cae1621f1221f9d4c6a9745dd7233f29de581b6", + "sha256:7169a34971e398dd58e87e173f97366fd88a3fa80852704530433eb224a8ca57", + "sha256:91c54d6bb9eeaaff965656c5ea6cbdcbf780bad8462ac99b30b451548194746f", + "sha256:aeef177647bb3fccfe09065481989d7dfc5ac59e9367d6a00a3481062cf651e4", + "sha256:cf8ae10559a78aee0409ede1e9d4fda03895433eeafe609dd9ed67e45f552db0", + "sha256:d51d0889d1c4d51c51a9822265c0494ea3e70a52bdd88358e0863daca46fa23a", + "sha256:de5ccd3500247f85fe4f9fad90f80a8bd397e4f110a4c33fabf95f07403e8372", + "sha256:e1d33589e32f482d0a7d1957bf473d43341115d40d33f578dad44432e47df7b7", + "sha256:e8d1939262aa6b36d0c51f50a50a43a04b9618d20db31e6c0192b1463067aeef", + "sha256:e918d51b1fda82a65fdf52d2f3914b2246481cc2a9cd10e223e6be6078916ff3" + ], + "index": "pypi", + "version": "==3.0.3" + }, + "mpmath": { + "hashes": [ + "sha256:fc17abe05fbab3382b61a123c398508183406fa132e0223874578e20946499f6" + ], + "version": "==1.1.0" + }, + "numpy": { + "hashes": [ + "sha256:1980f8d84548d74921685f68096911585fee393975f53797614b34d4f409b6da", + "sha256:22752cd809272671b273bb86df0f505f505a12368a3a5fc0aa811c7ece4dfd5c", + "sha256:23cc40313036cffd5d1873ef3ce2e949bdee0646c5d6f375bf7ee4f368db2511", + "sha256:2b0b118ff547fecabc247a2668f48f48b3b1f7d63676ebc5be7352a5fd9e85a5", + "sha256:3a0bd1edf64f6a911427b608a894111f9fcdb25284f724016f34a84c9a3a6ea9", + "sha256:3f25f6c7b0d000017e5ac55977a3999b0b1a74491eacb3c1aa716f0e01f6dcd1", + "sha256:4061c79ac2230594a7419151028e808239450e676c39e58302ad296232e3c2e8", + "sha256:560ceaa24f971ab37dede7ba030fc5d8fa173305d94365f814d9523ffd5d5916", + "sha256:62be044cd58da2a947b7e7b2252a10b42920df9520fc3d39f5c4c70d5460b8ba", + "sha256:6c692e3879dde0b67a9dc78f9bfb6f61c666b4562fd8619632d7043fb5b691b0", + "sha256:6f65e37b5a331df950ef6ff03bd4136b3c0bbcf44d4b8e99135d68a537711b5a", + "sha256:7a78cc4ddb253a55971115f8320a7ce28fd23a065fc33166d601f51760eecfa9", + "sha256:80a41edf64a3626e729a62df7dd278474fc1726836552b67a8c6396fd7e86760", + "sha256:893f4d75255f25a7b8516feb5766c6b63c54780323b9bd4bc51cdd7efc943c73", + "sha256:972ea92f9c1b54cc1c1a3d8508e326c0114aaf0f34996772a30f3f52b73b942f", + "sha256:9f1d4865436f794accdabadc57a8395bd3faa755449b4f65b88b7df65ae05f89", + "sha256:9f4cd7832b35e736b739be03b55875706c8c3e5fe334a06210f1a61e5c2c8ca5", + "sha256:adab43bf657488300d3aeeb8030d7f024fcc86e3a9b8848741ea2ea903e56610", + "sha256:bd2834d496ba9b1bdda3a6cf3de4dc0d4a0e7be306335940402ec95132ad063d", + "sha256:d20c0360940f30003a23c0adae2fe50a0a04f3e48dc05c298493b51fd6280197", + "sha256:d3b3ed87061d2314ff3659bb73896e622252da52558f2380f12c421fbdee3d89", + "sha256:dc235bf29a406dfda5790d01b998a1c01d7d37f449128c0b1b7d1c89a84fae8b", + "sha256:fb3c83554f39f48f3fa3123b9c24aecf681b1c289f9334f8215c1d3c8e2f6e5b" + ], + "index": "pypi", + "version": "==1.16.2" + }, + "pyparsing": { + "hashes": [ + "sha256:1873c03321fc118f4e9746baf201ff990ceb915f433f23b395f5580d1840cb2a", + "sha256:9b6323ef4ab914af344ba97510e966d64ba91055d6b9afa6b30799340e89cc03" + ], + "version": "==2.4.0" + }, + "pyqt5": { + "hashes": [ + "sha256:020e8771ad66429587d3db6cceb028be3435e2c65f417dcc8e3b3e644b6ab1d7", + "sha256:399fd9c3b3786210d9d2a86bc73255cc03c997b454480b7c0daf3e1a09e1ab58", + "sha256:c0565737d1d5fd40a7e14b6692bfec1f1e1f0be082ae70d1776c44561980e3c4", + "sha256:d4e88208dfd017636e4b1f407990d0c3b6cf47afed6be4f2fb6ca887ef513e4b" + ], + "index": "pypi", + "version": "==5.12.1" + }, + "pyqt5-sip": { + "hashes": [ + "sha256:0aed4266e52f4d11947439d9df8c57d4578e92258ad0d70645f9f69b627b35c7", + "sha256:15750a4a3718c68ba662eede7310612aedad70c9727497908bf2eeddd255d16f", + "sha256:2071265b7a603354b624b6cfae4ea032719e341799fb925961fd192f328fd203", + "sha256:31a59f76168df007b480b9382256c20f8898c642e1394df2990559f0f6389f66", + "sha256:39ce455bf5be5c1b20db10d840a536f70c3454be3baa7adca22f67ed90507853", + "sha256:5d8a4b8e75fa6c9f837ea92793f1aeacbf3929f1d11cdf36f9604c1876182aca", + "sha256:6bf9bacac776257390dd3f8a703d12c8228e4eb1aa0b191f1550965f44e4c23a", + "sha256:87a76841e07e6dae5df34dded427708b5ea8d605c017e0bb6e56a7354959271e", + "sha256:9f0ad4198f2657da9d59a393c266959362c517445cace842051e5ad564fa8fb0", + "sha256:cbb0f24e3bfa1403bf3a05f30066d4457509aa411d6c6823d02b0508b787b8ea", + "sha256:d38d379fa4f5c86a8ea92b8fab42748f37de093d21a25c34b2339eb2764ee5d5", + "sha256:f24413d34707525a0fe1448b9574089242b9baa7c604ae5714abef9ad1749839" + ], + "version": "==4.19.15" + }, + "python-dateutil": { + "hashes": [ + "sha256:7e6584c74aeed623791615e26efd690f29817a27c73085b78e4bad02493df2fb", + "sha256:c89805f6f4d64db21ed966fda138f8a5ed7a4fdbc1a8ee329ce1b74e3c74da9e" + ], + "version": "==2.8.0" + }, + "scipy": { + "hashes": [ + "sha256:014cb900c003b5ac81a53f2403294e8ecf37aedc315b59a6b9370dce0aa7627a", + "sha256:281a34da34a5e0de42d26aed692ab710141cad9d5d218b20643a9cb538ace976", + "sha256:588f9cc4bfab04c45fbd19c1354b5ade377a8124d6151d511c83730a9b6b2338", + "sha256:5a10661accd36b6e2e8855addcf3d675d6222006a15795420a39c040362def66", + "sha256:628f60be272512ca1123524969649a8cb5ae8b31cca349f7c6f8903daf9034d7", + "sha256:6dcc43a88e25b815c2dea1c6fac7339779fc988f5df8396e1de01610604a7c38", + "sha256:70e37cec0ac0fe95c85b74ca4e0620169590fd5d3f44765f3c3a532cedb0e5fd", + "sha256:7274735fb6fb5d67d3789ddec2cd53ed6362539b41aa6cc0d33a06c003aaa390", + "sha256:78e12972e144da47326958ac40c2bd1c1cca908edc8b01c26a36f9ffd3dce466", + "sha256:790cbd3c8d09f3a6d9c47c4558841e25bac34eb7a0864a9def8f26be0b8706af", + "sha256:79792c8fe8e9d06ebc50fe23266522c8c89f20aa94ac8e80472917ecdce1e5ba", + "sha256:865afedf35aaef6df6344bee0de391ee5e99d6e802950a237f9fb9b13e441f91", + "sha256:870fd401ec7b64a895cff8e206ee16569158db00254b2f7157b4c9a5db72c722", + "sha256:963815c226b29b0176d5e3d37fc9de46e2778ce4636a5a7af11a48122ef2577c", + "sha256:9726791484f08e394af0b59eb80489ad94d0a53bbb58ab1837dcad4d58489863", + "sha256:9de84a71bb7979aa8c089c4fb0ea0e2ed3917df3fb2a287a41aaea54bbad7f5d", + "sha256:b2c324ddc5d6dbd3f13680ad16a29425841876a84a1de23a984236d1afff4fa6", + "sha256:b86ae13c597fca087cb8c193870507c8916cefb21e52e1897da320b5a35075e5", + "sha256:ba0488d4dbba2af5bf9596b849873102d612e49a118c512d9d302ceafa36e01a", + "sha256:d78702af4102a3a4e23bb7372cec283e78f32f5573d92091aa6aaba870370fe1", + "sha256:def0e5d681dd3eb562b059d355ae8bebe27f5cc455ab7c2b6655586b63d3a8ea", + "sha256:e085d1babcb419bbe58e2e805ac61924dac4ca45a07c9fa081144739e500aa3c", + "sha256:e2cfcbab37c082a5087aba5ff00209999053260441caadd4f0e8f4c2d6b72088", + "sha256:e742f1f5dcaf222e8471c37ee3d1fd561568a16bb52e031c25674ff1cf9702d5", + "sha256:f06819b028b8ef9010281e74c59cb35483933583043091ed6b261bb1540f11cc", + "sha256:f15f2d60a11c306de7700ee9f65df7e9e463848dbea9c8051e293b704038da60", + "sha256:f31338ee269d201abe76083a990905473987371ff6f3fdb76a3f9073a361cf37", + "sha256:f6b88c8d302c3dac8dff7766955e38d670c82e0d79edfc7eae47d6bb2c186594" + ], + "index": "pypi", + "version": "==1.2.1" + }, + "six": { + "hashes": [ + "sha256:3350809f0555b11f552448330d0b52d5f24c91a322ea4a15ef22629740f3761c", + "sha256:d16a0141ec1a18405cd4ce8b4613101da75da0e9a7aec5bdd4fa804d0e0eba73" + ], + "version": "==1.12.0" + }, + "sympy": { + "hashes": [ + "sha256:71a11e5686ae7ab6cb8feb5bd2651ef4482f8fd43a7c27e645a165e4353b23e1", + "sha256:f9b00ec76151c98470e84f1da2d7d03633180b71fb318428ddccce1c867d3eaa" + ], + "index": "pypi", + "version": "==1.4" + } + }, + "develop": {} +} diff --git a/saved_u.pickle b/saved_u.pickle new file mode 100644 index 0000000..06ca3fd Binary files /dev/null and b/saved_u.pickle differ diff --git a/src/__pycache__/solver.cpython-37.pyc b/src/__pycache__/solver.cpython-37.pyc deleted file mode 100644 index 15bef4d..0000000 Binary files a/src/__pycache__/solver.cpython-37.pyc and /dev/null differ diff --git a/src/app.py b/src/app.py index a073d5e..f0114cc 100644 --- a/src/app.py +++ b/src/app.py @@ -1,11 +1,14 @@ import sys from PyQt5.QtWidgets import QApplication, QGridLayout, QGroupBox, QLabel, QLineEdit, QMessageBox, QPushButton, \ - QSizePolicy, QVBoxLayout, QWidget + QSizePolicy, QVBoxLayout, QWidget, QCheckBox from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas from matplotlib.figure import Figure -import solver +from solvers.crank_nicolson import CrankNicolsonSolver +from solvers.explicit import ExplicitSolver +from solvers.implicit import ImplicitSolver +from solvers.theoretical import TheoreticalSolver class App(QWidget): @@ -13,7 +16,16 @@ class App(QWidget): def __init__(self): super().__init__() - self.solver = solver.Solver() + self.theoretical_solver = TheoreticalSolver() + self.explicit_solver = ExplicitSolver() + self.implicit_solver = ImplicitSolver() + self.crank_nicolson_solver = CrankNicolsonSolver() + self.solvers = [ + self.theoretical_solver, + self.explicit_solver, + self.implicit_solver, + self.crank_nicolson_solver, + ] self.grid = QGridLayout() self.setLayout(self.grid) @@ -29,34 +41,57 @@ class App(QWidget): param_vbox.addWidget(QLabel('Длина цилиндра l')) self.l_field = QLineEdit('20') - self.l_field.setStyleSheet("color: black;") + # self.l_field.setStyleSheet("color: black;") param_vbox.addWidget(self.l_field) param_vbox.addWidget(QLabel('Коэффициент диффузии D')) self.D_field = QLineEdit('0.6') - self.D_field.setStyleSheet("color: black;") + # self.D_field.setStyleSheet("color: black;") param_vbox.addWidget(self.D_field) param_vbox.addWidget(QLabel('Мембранный коэффициент диффузии H')) self.H_field = QLineEdit('4') - self.H_field.setStyleSheet("color: black;") + # self.H_field.setStyleSheet("color: black;") param_vbox.addWidget(self.H_field) param_vbox.addWidget(QLabel('Момент времени t')) self.t_field = QLineEdit('0') - self.t_field.setStyleSheet("color: black;") + # self.t_field.setStyleSheet("color: black;") param_vbox.addWidget(self.t_field) param_vbox.addWidget(QLabel('Точность поиска собственных значений')) self.eps_field = QLineEdit('1e-9') - self.eps_field.setStyleSheet("color: black;") + # self.eps_field.setStyleSheet("color: black;") param_vbox.addWidget(self.eps_field) param_vbox.addWidget(QLabel('Количество членов ряда в суммировании')) self.sum_count_field = QLineEdit('100') - self.sum_count_field.setStyleSheet("color: black;") + # self.sum_count_field.setStyleSheet("color: black;") param_vbox.addWidget(self.sum_count_field) + param_vbox.addWidget(QLabel('Размерность сетки по времени K')) + self.K_field = QLineEdit('100') + # self.K_field.setStyleSheet("color: black;") + param_vbox.addWidget(self.K_field) + + param_vbox.addWidget(QLabel('Размерность сетки по координате I')) + self.I_field = QLineEdit('100') + # self.I_field.setStyleSheet("color: black;") + param_vbox.addWidget(self.I_field) + + self.theoretical_checkbox = QCheckBox('Аналитическое решение') + self.theoretical_checkbox.setChecked(True) + param_vbox.addWidget(self.theoretical_checkbox) + self.explicit_checkbox = QCheckBox('Явная схема') + self.explicit_checkbox.setChecked(True) + param_vbox.addWidget(self.explicit_checkbox) + self.implicit_checkbox = QCheckBox('Неявная схема') + self.implicit_checkbox.setChecked(True) + param_vbox.addWidget(self.implicit_checkbox) + self.crank_nicolson_checkbox = QCheckBox('Схема Кранка-Николсон') + self.crank_nicolson_checkbox.setChecked(True) + param_vbox.addWidget(self.crank_nicolson_checkbox) + calc_button = QPushButton('Рассчитать') calc_button.pressed.connect(self.calc) param_vbox.addWidget(calc_button) @@ -121,22 +156,64 @@ class App(QWidget): self.show_error(f'Некорректный формат числа: {sum_count_str}') return + try: + K_str = self.K_field.text() + K_val = int(K_str) + except Exception as e: + self.show_error(f'Некорректный формат числа: {K_str}') + return + + try: + I_str = self.I_field.text() + I_val = int(I_str) + except Exception as e: + self.show_error(f'Некорректный формат числа: {I_str}') + return + if not ( - 1e-3 <= l_val <= 1e+3 and - 1e-3 <= D_val <= 1e+3 and - 1e-3 <= H_val <= 1e+3 and - 0 <= t_val <= 1e+3 and - 1e-20 <= eps_val <= 1e+3 and - 1 <= sum_count_val <= 1e+3 + 1e-3 <= l_val <= 1e+3 and + 1e-3 <= D_val <= 1e+3 and + 1e-3 <= H_val <= 1e+3 and + 0 <= t_val <= 1e+3 and + 1e-20 <= eps_val <= 1e+3 and + 1 <= sum_count_val <= 1e+3 and + 1 <= K_val <= 1e+4 and + 1 <= I_val <= 1e+4 ): self.show_error(f'Некорректные входные данные') return - self.solver.set_params(l_val, D_val, H_val, t_val, eps_val, sum_count_val) + for solver in self.solvers: + solver.set_params( + l=l_val, + D=D_val, + H=H_val, + t=t_val, + eps=eps_val, + sum_count=sum_count_val, + K=K_val, + I=I_val, + ) - plot_x, plot_y = self.solver.calc() + self.plot.prep() - self.plot.plot(plot_x, plot_y) + if self.theoretical_checkbox.isChecked(): + plot_x, plot_y = self.theoretical_solver.calc() + self.plot.plot(plot_x, plot_y, label='Аналитическое решение', linestyle='-') + + if self.explicit_checkbox.isChecked(): + plot_x, plot_y = self.explicit_solver.calc() + self.plot.plot(plot_x, plot_y, label='Явная схема', linestyle='--') + + if self.implicit_checkbox.isChecked(): + plot_x, plot_y = self.implicit_solver.calc() + self.plot.plot(plot_x, plot_y, label='Неявная схема', linestyle='-.') + + if self.crank_nicolson_checkbox.isChecked(): + plot_x, plot_y = self.crank_nicolson_solver.calc() + self.plot.plot(plot_x, plot_y, label='Схема Кранка-Николсон', linestyle=':') + + self.plot.show() class PlotCanvas(FigureCanvas): @@ -150,16 +227,21 @@ class PlotCanvas(FigureCanvas): FigureCanvas.setSizePolicy(self, QSizePolicy.Expanding, QSizePolicy.Expanding) FigureCanvas.updateGeometry(self) - def plot(self, plot_x, plot_y): + def prep(self): self.axes.clear() - self.axes.plot(plot_x, plot_y) + def plot(self, plot_x, plot_y, **kwargs): + self.axes.plot(plot_x, plot_y, **kwargs) + + def show(self): self.axes.set_xlabel('$x, \\mathrm{м}$') self.axes.set_ylabel('$u$') self.axes.set_title('Поле концентрации вещества в цилиндре, $\\frac{\\mathrm{кг}}{\\mathrm{м}^3}$') self.axes.set_ylim(0, None) + self.axes.legend() + self.axes.grid(True) self.draw() diff --git a/src/solvers/__init__.py b/src/solvers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/solvers/base.py b/src/solvers/base.py new file mode 100644 index 0000000..67caeef --- /dev/null +++ b/src/solvers/base.py @@ -0,0 +1,19 @@ +from typing import Tuple +import numpy as np + + +class Solver: + def __init__(self): + self.l = None + self.D = None + self.H = None + self.t = None + + def set_params(self, **kwargs): + self.l = kwargs['l'] + self.D = kwargs['D'] + self.H = kwargs['H'] + self.t = kwargs['t'] + + def calc(self) -> Tuple[np.ndarray, np.ndarray]: + raise NotImplementedError() diff --git a/src/solvers/crank_nicolson.py b/src/solvers/crank_nicolson.py new file mode 100644 index 0000000..0bce039 --- /dev/null +++ b/src/solvers/crank_nicolson.py @@ -0,0 +1,76 @@ +from typing import Tuple + +import numpy as np +import scipy.sparse.linalg + +from solvers.base import Solver +from util import Grid + + +class CrankNicolsonSolver(Solver): + def __init__(self): + super().__init__() + + self.K = None + self.I = None + + self.params = None + + def set_params(self, **kwargs): + super().set_params(**kwargs) + + self.K = kwargs['K'] + self.I = kwargs['I'] + + self.params = kwargs + + def calc(self) -> Tuple[np.ndarray, np.ndarray]: + grid = Grid(**self.params) + + vals = solve_crank_nicolson(grid) + + return grid.xs, vals[-1, :] + + +# @cachier(stale_after=datetime.timedelta(hours=1)) +def solve_crank_nicolson(grid: Grid) -> np.ndarray: + ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float) + + ug[0, :] = np.sin(np.pi * grid.xs / grid.l) ** 2 * 4 / 10 + + for k in range(1, grid.K + 1): + diags = { + -1: np.zeros(grid.I + 1, dtype=float), + 0: np.zeros(grid.I + 1, dtype=float), + +1: np.zeros(grid.I + 1, dtype=float), + } + b = np.zeros(grid.I + 1, dtype=float) + + # left + diags[0][0] = grid.H * grid.h_x + 1 + diags[+1][0] = -1 + + # middle + for i in range(1, grid.I): + diags[-1][i] = grid.gamma + diags[0][i] = - 2 * grid.gamma - 2 + diags[+1][i] = grid.gamma + b[i] = -( + ug[k - 1][i - 1] * grid.gamma + + ug[k - 1][i] * (-2 * grid.gamma + 2) + + ug[k - 1][i + 1] * grid.gamma + ) + + # right + diags[0][grid.I] = grid.H * grid.h_x + 1 + diags[-1][grid.I] = -1 + + A = scipy.sparse.diags([ + diags[-1][1:], + diags[0], + diags[+1][:-1], + ], [-1, 0, 1], format='csc') + + ug[k, :] = scipy.sparse.linalg.spsolve(A, b) + + return ug diff --git a/src/solvers/explicit.py b/src/solvers/explicit.py new file mode 100644 index 0000000..ed5150e --- /dev/null +++ b/src/solvers/explicit.py @@ -0,0 +1,57 @@ +from typing import Optional, Tuple + +import numpy as np + +from solvers.base import Solver +from util import Grid + + +class ExplicitSolver(Solver): + def __init__(self): + super().__init__() + + self.K = None + self.I = None + + self.params = None + + def set_params(self, **kwargs): + super().set_params(**kwargs) + + self.K = kwargs['K'] + self.I = kwargs['I'] + + self.params = kwargs + + def calc(self) -> Tuple[np.ndarray, np.ndarray]: + grid = Grid(**self.params) + + vals = solve_explicit(grid) + + return grid.xs, vals[-1, :] + + +# @cachier(stale_after=datetime.timedelta(hours=1)) +def solve_explicit(grid: Grid) -> Optional[np.ndarray]: + # if grid.gamma > 1 / 2: + # return None + + ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float) + + ug[0, :] = np.sin(np.pi * grid.xs / grid.l) ** 2 * 4 / 10 + + for k in range(1, grid.K + 1): + # middle + ug[k, 1: grid.I] = ( + grid.gamma * ug[k - 1, 0: grid.I - 1] + + (-2 * grid.gamma + 1) * ug[k - 1, 1: grid.I] + + grid.gamma * ug[k - 1, 2: grid.I + 1] + ) + + # left + ug[k, 0] = ug[k, 1] / (grid.H * grid.h_x + 1) + + # right + ug[k, grid.I] = ug[k, grid.I - 1] / (grid.H * grid.h_x + 1) + + return ug diff --git a/src/solvers/implicit.py b/src/solvers/implicit.py new file mode 100644 index 0000000..62b4859 --- /dev/null +++ b/src/solvers/implicit.py @@ -0,0 +1,72 @@ +from typing import Optional, Tuple + +import numpy as np +import scipy.sparse.linalg + +from solvers.base import Solver +from util import Grid + + +class ImplicitSolver(Solver): + def __init__(self): + super().__init__() + + self.K = None + self.I = None + + self.params = None + + def set_params(self, **kwargs): + super().set_params(**kwargs) + + self.K = kwargs['K'] + self.I = kwargs['I'] + + self.params = kwargs + + def calc(self) -> Tuple[np.ndarray, np.ndarray]: + grid = Grid(**self.params) + + vals = solve_implicit(grid) + + return grid.xs, vals[-1, :] + + +# @cachier(stale_after=datetime.timedelta(hours=1)) +def solve_implicit(grid) -> np.ndarray: + ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float) + + ug[0, :] = np.sin(np.pi * grid.xs / grid.l) ** 2 * 4 / 10 + + for k in range(1, grid.K + 1): + diags = { + -1: np.zeros(grid.I + 1, dtype=float), + 0: np.zeros(grid.I + 1, dtype=float), + +1: np.zeros(grid.I + 1, dtype=float), + } + b = np.zeros(grid.I + 1, dtype=float) + + # left + diags[0][0] = grid.H * grid.h_x + 1 + diags[+1][0] = -1 + + # middle + for i in range(1, grid.I): + diags[-1][i] = grid.gamma + diags[0][i] = - 2 * grid.gamma - 1 + diags[+1][i] = grid.gamma + b[i] = -ug[k - 1, i] + + # right + diags[-1][grid.I] = -1 + diags[0][grid.I] = grid.H * grid.h_x + 1 + + A = scipy.sparse.diags([ + diags[-1][1:], + diags[0], + diags[+1][:-1], + ], [-1, 0, 1], format='csc') + + ug[k, :] = scipy.sparse.linalg.spsolve(A, b) + + return ug diff --git a/src/solver.py b/src/solvers/theoretical.py similarity index 88% rename from src/solver.py rename to src/solvers/theoretical.py index 457fa31..d90a98a 100644 --- a/src/solver.py +++ b/src/solvers/theoretical.py @@ -3,13 +3,12 @@ import pickle import numpy as np +from solvers.base import Solver -class Solver: + +class TheoreticalSolver(Solver): def __init__(self): - self.l = None - self.D = None - self.H = None - self.t = None + super().__init__() self.eps = 1e-9 self.sum_count = 100 @@ -34,13 +33,11 @@ class Solver: else: l = m - def set_params(self, l, D, H, t, eps, sum_count): - self.l = l - self.D = D - self.H = H - self.t = t - self.eps = eps - self.sum_count = sum_count + def set_params(self, **kwargs): + super().set_params(**kwargs) + + self.eps = kwargs['eps'] + self.sum_count = kwargs['sum_count'] self.alphas = np.array([self.find_root(i) for i in range(self.sum_count)]) @@ -109,6 +106,6 @@ def u(x, t): if __name__ == '__main__': - solver = Solver() - solver.set_params(20, 0.6, 4, 10) - print(solver.calc(10)) + solver = TheoreticalSolver() + # solver.set_params(20, 0.6, 4, 10) + # print(solver.calc(10)) diff --git a/src/util.py b/src/util.py new file mode 100644 index 0000000..3d34382 --- /dev/null +++ b/src/util.py @@ -0,0 +1,54 @@ +from dataclasses import dataclass, field + +import numpy as np + + +@dataclass +class Grid: + I: int + K: int + + xs: np.ndarray = field(repr=False, compare=False) + ts: np.ndarray = field(repr=False, compare=False) + + h_x: float = field(repr=False, compare=False) + h_t: float = field(repr=False, compare=False) + + gamma: float = field(repr=False, compare=False) + + l: float + t: float + D: float + H: float + + def __init__(self, I: int, K: int, l: float, t: float, D: float, H: float, **kwargs): + self.I = I + self.K = K + + self.xs = np.linspace(0, l, I + 1) + self.ts = np.linspace(0, t, K + 1) + + self.h_x = l / I + self.h_t = t / K + + self.gamma = D * self.h_t / (self.h_x ** 2) + + self.l = l + self.t = t + self.D = D + self.H = H + + def __hash__(self): + return hash((self.I, self.K)) + + +def mse(ug1: np.ndarray, ug2: np.ndarray) -> float: + return np.mean((ug1 - ug2) ** 2) + + +def mae(ug1: np.ndarray, ug2: np.ndarray) -> float: + return np.max(np.abs(ug1 - ug2)) + + +def mae_last(ug1: np.ndarray, ug2: np.ndarray) -> float: + return np.max(np.abs(ug1[-1, :] - ug2[-1, :]))