Блог

Главная - Вопросы по программированию - numpy cumsum с граничными условиями

numpy cumsum с граничными условиями

#python #numpy #cumsum

#python #numpy #cumsum

Вопрос:

Мне нужно выполнить две интеграции, чтобы получить совокупный перенос воды в области на основе профилей плотности и ширины между профилями.

Во-первых, скорость вычисляется как;

Обратите внимание, что V может быть ненулевым на поверхности и должно быть нулевым на глубине H.

Во-вторых, транспорт вычисляется как;

В этой интеграции T должно быть равно нулю на поверхности и внизу.

На данный момент я закодировал приведенные выше уравнения, используя numpy.ma.cumsum

 # Compute velocity
V = constant*numpy.ma.cumsum((density/width*dz)[::-1],axis=0)[::-1]
# Compute transport in Sverdrup (10**6 m/s^2)
T = width * numpy.ma.cumsum((V[::-1]-1/H*numpy.ma.cumsum((V*dz)[::-1],axis=0))*dz[::-1],axis=0)[::-1] / 10**6

Эти уравнения дают мне следующие результаты;
профиль плотности, скорость и транспортировка

Если у вас есть идея о том, как заставить транспорт T стремиться к нулю на поверхности и на глубине H, пожалуйста, дайте мне знать.

Вот мои наборы данных:

 constant = -100

density =      [0.08162372168703769, -0.04471999854399655,
               -0.3110550684352802, -0.4783609962985338,
               -0.2606431434960541, -0.1734227201536669,
               -0.20306011069737906, -0.22652204210817217,
               -0.24526932281150948, -0.2899378668189456,
               -0.2968024166812029, -0.30554869432489795,
               -0.3145306210531089, 0.04810820373199931,
               0.08374671599813155, 0.11754102251029508,
               0.14828751252093753, 0.17433953329941687,
               0.19374034901125015, 0.20420143703017857,
               0.2026197673765182, 0.1848936664282519,
               0.131398194555004, 0.3103135616245254,
               0.3938907065328294, 0.4770966560975012,
               0.5276251000107095, 0.48519935265176173,
               0.41467987590567645, 0.34373320113900263,
               0.12321804603743658, 0.05310501021676828,
               0.006689604194434651, -0.02221187988652673,
               -0.03928717776148005, -0.06361049772976912,
               -0.06657255480854474, -0.07648006375347904,
               -0.0671293538559894, -0.04341857793974668,
               -0.04089445488693855, -0.025967806515836855,
               0.010263609448429634, 0.016120765154255423,
               0.016567011411552812, 0.02229294316521191,
               0.02531710921425656, 0.027188241842850402,
               0.027440846175068145, --]

width =       [5370872.025103143, 5370872.025103143,
               5370872.025103143, 5370872.025103143,
               5286952.149710908, 5203032.274318671,
               5203032.274318671, 5203032.274318671,
               5203032.274318671, 5035192.523534197,
               5035192.523534197, 5035192.523534197,
               5035192.523534197, 4867352.772749724,
               4867352.772749724, 4867352.772749724,
               4867352.772749724, 4867352.772749724,
               4867352.772749724, 4867352.772749724,
               4867352.772749724, 4867352.772749724,
               4867352.772749724, 4783432.897357487,
               4783432.897357487, 4783432.897357487,
               4783432.897357487, 4699513.02196525, 4699513.02196525,
               4699513.02196525, 4699513.02196525, 4699513.02196525,
               4699513.02196525, 4699513.02196525, 4699513.02196525,
               4699513.02196525, 4615593.146573015, 4615593.146573015,
               4615593.146573015, 4615593.146573015,
               4531673.271180778, 4531673.271180778,
               4447753.395788541, 4447753.395788541,
               4363833.520396304, 4279913.645004068,
               4195993.769611831, 4028154.018827358,
               3776394.3926506485, --]

dz =          [-5.03355, -10.0671, -10.118699999999999,
               -10.139100000000006, -10.217899999999986,
               -10.276900000000012, -10.40849999999999,
               -10.541100000000029, -10.774099999999947,
               -11.046300000000059, -11.472899999999953,
               -12.010500000000036, -12.801899999999975,
               -13.838099999999997, -15.307100000000048,
               -17.25869999999992, -19.95050000000009,
               -23.514299999999935, -28.290900000000022,
               -34.48270000000002, -42.456899999999905,
               -52.35730000000012, -64.36129999999991, -78.2561,
               -93.72070000000008, -109.98889999999983,
               -126.30150000000026, -141.72169999999983,
               -155.66089999999986, -167.62050000000045,
               -177.5652999999993, -185.49530000000095,
               -191.7308999999991, -196.46390000000042,
               -200.07110000000011, -202.71990000000005,
               -204.71409999999923, -206.14090000000124,
               -207.21909999999843, -207.96890000000167,
               -208.54709999999795, -208.93390000000272,
               -209.24569999999676, -209.43890000000374,
               -209.61189999999624, -209.70250000000306,
               -209.80249999999796, -209.84110000000146,
               -209.90109999999913, --]

H = 5185.11075

integral of dx = width

Комментарии:

1. Как x входит в уравнение для T(z) ? Ничто в этом уравнении не параметризуется x — это просто интеграл от константы?

2. Уравнение для T(z) не показывает ничего, что указывает на это T(0) = T(-H) = 0 . Есть ли какие-то дополнительные отношения, которые здесь не показаны?

3. Единственное, что я вижу в этих уравнениях, которыми (потенциально) можно манипулировать, — это значение constant . Я не вижу, как это может обеспечить соблюдение ваших граничных условий, и я не вижу ничего другого, что можно было бы скорректировать. Похоже, это скорее математическая проблема, чем проблема numpy

4. Верхний предел и переменная интегрирования должны отличаться в обоих интегралах. Не могли бы вы пометить свои графики так, чтобы они соответствовали используемым вами переменным?

5. Почему вы используете np.ma , когда у вас, похоже, нет замаскированного массива? Или ты? Не могли бы вы опубликовать MCVE?