This commit is contained in:
2019-04-12 13:32:38 +04:00
parent 5c01c4d176
commit 9a2360d13b
13 changed files with 603 additions and 35 deletions

1
.gitignore vendored Normal file
View File

@@ -0,0 +1 @@
.venv

View File

@@ -4,6 +4,11 @@ verify_ssl = true
name = "pypi"
[packages]
"pyqt5" = "*"
sympy = "*"
numpy = "*"
scipy = "*"
matplotlib = "*"
[dev-packages]

205
Pipfile.lock generated Normal file
View File

@@ -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": {}
}

BIN
saved_u.pickle Normal file

Binary file not shown.

View File

@@ -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()

0
src/solvers/__init__.py Normal file
View File

19
src/solvers/base.py Normal file
View File

@@ -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()

View File

@@ -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

57
src/solvers/explicit.py Normal file
View File

@@ -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

72
src/solvers/implicit.py Normal file
View File

@@ -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

View File

@@ -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))

54
src/util.py Normal file
View File

@@ -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, :]))