Здравствуйте!
Я вот решаю квазилинейное уравнение теплопроводности, а именно dT/dt = k d/dx ( T^s *dT/dx ). Прошу прощения, не знаю как здесь нормально вставить формулу...
В общем, проблема в том, что обычная разностная схема не работает. Вот такая:

где Т2 - температура в данный момент времени, oldT2 - температура в предыдущий момент времени, dt - промежуток по времени, dh - промежуток по координате. А именно, температура вовсе не распространяется.
Гугл показывал мне Самарского, вроде там что-то на эту тему есть, но я не нашла. То есть поняла, что существует такой интегро-итерационый метод (если я правильно запомнила его название) но он меня стабильно приводит к системе нелинейных уравнений, которую я решать вообще не знаю как.
Посоветуйте что-нибудь, пожалуйста!

Что такое "обычный разностный метод"? Подсчет тепловой волны с помощью т.н. "наилучшей аппроксимации"? Он и не должен работать, это нормально. Воспользуйтесь аппроксимацией потока в виде J_n= \frac{T^s_{n-1/2}+T^s_{n+1/2}}{2} \cdot \frac{T_{n+1/2}-T_{n-1/2}}{(h_{n-1/2}+h_{n+1/2})/2}. Для степенного коэффициента теплопроводности эта аппроксимация точна.

Фантом, спасибо за подсказку!
А можно спросить пару вопросов, например, что такое J_n ? Температура в следующий момент времени?
И потом, T_{n+1/2}, это имеется в виду, индекс по координате? И температуру там считать как среднее от температуры в двух узлах?

Поток в узле n. Т.е. то, что находится внутри "внешней" производной.


Да. В таких задачах обычно используют "полуцелые" сетки, так удобнее.

Понятно, спасибо! Попробую.
И еще вот: (h_{n-1/2}+h_{n+1/2})/2 - это зачем? Это если расстояние между узлами не одинаковое? Я правильно поняла, что h - это расстояние между узлами по координате?

Это сообщение отредактировал(а) Kafari - 3.4.2011, 22:54

Да. Поскольку я не знаю, используете Вы равномерную сетку или нет, на всякий случай написал более общий вариант.

Наверное, я жестоко туплю, но у меня все равно ничего не работает. Код выглядит так вот:


Помогите найти ошибку, пожалуйста!

на вскидку (остальное не смотрел): 1/2 получается 0, т.к. деление будет целочисленное
замени либо на 0.5 либо 1.0/2.0 как больше нравится

вообще все деления должны приводиться к double


Это сообщение отредактировал(а) borisbn - 4.4.2011, 17:07

Это, что, Ваша первая программа на C/С++?  Результат вычисления 1/2 - это нуль. Соответственно, J1 и J2 всегда окажутся нулевыми.

Заодно уж уберите двукратное вычисление pow(Ts2,s) и создание/удаление переменных в каждой ветке цикла. На численный результат это не повлияет, но вычислительные программы с таким "качеством" оптимизации писать неприлично.

Спасибо за подсказку! Исправила.


Вторая. До этого был матлаб, поэтому у меня еще бывают всякие глюки.. В любом случае спасибо, буду знать.
Теперь код выглядит так:

Теперь правильно?

В принципе, да, хотя на самом деле использование "полуцелой" сетки обычно предполагает и хранение данных с ней же. Температуру тут проще считать на полуцелых узлах, а потоки - на целых. Кроме этого, от двукратного вычисления pow(Ts2,s) я бы все-таки избавился - возведение в степень штука дорогая. Ну и делать постоянно приведение типов для dh тоже как-то странно - не проще ли сразу объявить ее как double?

Фантом, соптимизировать быстрее, чем Intel'овский компилятор или чем Intel Performance Primitives всё равно не получится. Перед тем как начинать гнаться за скоростью нужно полностью отладить алгоритм/логику/математику.
Мало того, создание переменной на стеке ("в каждой ветке цикла") быстрее, чем выделение в куче или даже в статике.
(С) Александреску - преждевременная оптимизация есть источник зла.
Думайте только о том, чтобы программа работала правильно.

Как-то патетически получилось... 

Хм... Ну что же, попробуем. Я приписал к этому куску какую-то обвязку (чтобы было что компилировать) и погонял в двух видах. В одном цикл был исходным, в другом - возведение в степень делалось один раз (из двух) и записывалось в промежуточную переменную, после чего результат два раза использовался. И то, и другое компилировалось ICC с ключом -fast. Прирост производительности, конечно, получился смешным (примерно 3%), но он систематический.  


Возможно, но, во-первых, после такой отладки "логики" найти ее в программе будет затруднительно, и, во-вторых, это действие попросту уменьшает вероятность случайной ошибки.


Да, тут согласен (хотя, вообще говоря, это до некоторой степени "компиляторозависимо"). Правда, там как минимум две переменных из пяти вообще не нужны (а при аккуратном расписывании логики удастся убрать и больше двух).


Если бы еще Александреску занимался вычислительным программированием...   

А это как? Сделать еще одну переменную?


Ясно, сделано.

Вообще я понимаю, что промежуточные переменные не нужны, и можно было бы сразу написать длинную строчку. Возможно, так и сделаю, но пока так понятнее.

Вот, а теперь у меня проблема с визуализацией... Делается у меня это в Qt, там рисуется полосочка (пока одномерный случай), разделенная на узкие прямоугольники, закрашиваемые определенным цветом. Крайняя справа и крайняя слева фиксированы - там температура постоянна. Должен получиться градиент, но он не получается. Собственно, вот как оно выглядит, если это понятнее, чем на словах объяснять

Если что, dataAvailable - это такая переменная, чтобы отрисовка включалась именно тогда, когда надо, а не просто когда вдруг update вызовется, мало ли...

Но то, что рисуется, слабо похоже на правду. То есть правильно получается только в одном случае - когда s (степень) равна нулю. В остальном я наблюдаю схлопывающиеся белые полосы. Вроде такого вот

Здесь нижняя часть окошка, вверху распространение нормального тепла...
И вот я думаю, то ли температура там становится отрицательная, то ли слишком большая наоборот, то ли у меня с визаулизатором проблема. Насколько я поняла из намеков компилятора, он пишет " "QColor::setAlpha": invalid value -2147483648 ". Что бы это могло значить?

Неправильно понимаешь. Промежуточные переменные нужны. Так проще отлаживаться. Так понятней. 

Ключевое слово - понятней. Код в разы чаще читается, чем компилируется...

Теперь к отрисовке:
1. Ты уверена, что у тебя получаются "правильные" данные ? Распечатай их в qDebug и посмотри.
2. "Насколько я поняла из намеков компилятора". Правильнее сказать не из компилятора, а из отладочного вывода Qt (компилятор уже отработал. программы выполняется). Так, вернёмся к нашим... 
     Посмотри (выведи в qDebug) чему равно (data->oldT2[x])*225/data->T0; если оно больше 255, то ессно ругается, что alpha < 0
3. Насколько я понял из алгоритма, ты устанавливаешь прозрачность чёрного цвета, в зависимости от значения функции. Может просто попробовать заливать цветом ?


В общем, попробуй это и покажи, что выводится в Output


Это сообщение отредактировал(а) borisbn - 5.4.2011, 08:40