Система громоздких ОДУ. Зашел в тупик, нужен пинок.

[uni]

Может в Mathcad есть еще какое-нибудь чудесные средства для поиска "мест деления на ноль".

Если их нет, то их можно придумать:

[uni]

Всякие такие вещи в книжках вряд ли найти можно:



Тут простой способ сокращения записи. Люблю, когда всё компактно.

[uni]

Вот что можно сделать таким подходом:



Профилирование кода документа Mathcad (pdf)

Исходник: MathcadReport.xmcdz

[uni]

Как это работает:

[Maxakch]

Вы не совсем разглядели мощь метода

я совсем не разглядел мощь этого метода, теперь со всем вроде разобрался. Я и не догадывался, что маткад так умеет. Ошибку пока не нашел. Все та же ошибка с слишком большим числом, теперь правда понятно, что она происходит при вызове функции f1, однако, в самой функции f1 деления на ноль не происходит. Значит нужно проверять те функции которые вызывает f1. Дома ноут очень слабый, но все равно попробую.
P.S. еще не разобрался с тем, как вы сделали вывод в pdf-ный документ. У меня окно трассировки обрезает, начало трассировки, т.е. выводится N-ное количество последних строк. Этого вполне достаточно для отладки, просто любопытно, как вы вывели весь документ.

[uni]

Код: Выделить всё

Mathcad 15.0 (M010 [MC15_M010_20110622])

Брал на рутрекере. Никаких ограничений по размеру отладочной информации у меня нет, наоборот, я даже обрезал лишнюю начальную внутреннюю отладку. Mathcad, перед тем как решать ДУ функцией rkfixed(), осуществляет какие-то внутренние манипуляции с вектором D(t,m), что я вижу. Только потом происходят непосредственно вычисления.

Прежде чем применять тяжёлую артиллерию в виде других решателей нужно по крайней мере получить хотя бы несколько точек стандартным фиксированным Рунге-Кутта. Поэтому я ограничил диапазон и количество точек до тривиальных значений, упрощая задачу.

Кроме того, я убрал кое-что в коде функции j2(), где под корнем, в отличие от j1(), был вызов функции Ф(). Именно на неё ругался решатель диффура, но я не знаю конечно правильно ли это. После такой правки я смог получить всё дерево вызовов для одного цикла работы решателя диффура. Но появилась проблема с комплексным результатом, от чего решатель стал выдавать новую ошибку: "результат должен быть реальным". Мнимая часть появилась как раз в том месте, где я исправлял. Это наглядно видно по отладочному листингу в pdf - там у некоторых результатов "i" появляется, что в конце концов приводит к комплексному результату.

Вариантов ошибок в таких ситуациях может быть несколько:

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

- комплексный результат (обычно это не ошибка) - очевидно, что это отрицательные значения под квадратным корнем где-то. Можно все такие места "пометить" при помощи trace(), выдавая сообщение в лог, если подкоренное выражение стало отрицательным. В принципе это и так видно по отладочному журналу в дереве вызовов - там результат функции имеет множитель i;

- ещё вариант - это глубокая рекурсия или зацикливание кода. Mathcad имеет ограничения по глубине вложенности и выдаст предупреждение, если "опознает" такую ситуацию.

Вот в общем-то и всё. Эти варианты ошибок проявляются только во время интерпретации кода. Остальные ошибки подсвечиваются красным при проверке целостности сборки документа, т.е. должны присутствовать все элементы, вызов которых имеется.

А pdf я делал в MS Word, просто скопировав код из окна отладочного журнала и напечатав его на принтере pdf.

Не забываем также, что кроме trace(), существуют ещё операторы: error(), on error(), return() и bp(). Последний делает то же, что и trace(), но приостанавливает интерпретацию, пока пользователь не нажмёт нужную клавишу на панели Отладки (есть такая панель).

П.С. Очень редко кто использует вывод отладочной информации таким образом. Такой метод использования отладочных возможностей Mathcad нужно брать на вооружение. Вряд ли даже на форуме PTC где-то можно увидеть такое.

[uni]

Покажите документ с исправленными формулами и я укажу место, где возникает переполнение.

[uni]

Я выше ошибся, не bp(), а pause(). bp() - так называлась моя аналогичная отладочная функция для моего отладчика.

[Maxakch]

Вообщем, семейные дела не дали мне поработать на выходных. Так, что опишу тут, что у меня получилось сегодня.

Кроме того, я убрал кое-что в коде функции j2(), где под корнем, в отличие от j1(), был вызов функции Ф(). Именно на неё ругался решатель диффура, но я не знаю конечно правильно ли это.

Собственно, я вернул вызов функции Ф()под этим корнем. И появилась ошибка с "слишком большим числом" в f1. Эту проблему я решил: Почему то Mathcad действительно не нравится гладкая Ф(), я заменил ее на функцию Ф(x)=if(x<0,0,1). Видимо в Ф(x)=1/(e^(-50*x)+1) происходит ошибка из-за экспоненты в степени.
Вообщем, теперь все считается (кроме Ф(x) поправил кое-что в коэффициентах). Но результаты мне не нравятся. Они мне кажутся не физичными, особенно в плане температуры, которая уменьшается. Ведь в моей задаче в нижнее помещение поступает разогретый пар, который с ростом концентрации может конденсировать. Т.е. "количество энергии" должно увеличивать в помещениях, что должно вести к увеличению температуры. ..

Да, и метод rkfixed и Radau дают разные решения, возможно нужно увеличить количество шагов, что бы они показывали одно и то же.

Так же попробовал запихнуть свою систему в Maple.Кстати, действительно удобнее вводить уравнения. И я так понял в Maple не обязательно представлять систему ОДУ в виде df(m)/dt=f(m), как это нужно делать в Mathcad. А это действительно удобнее, не нужно выписывать эти громоздкие f1...f6. Но maple выдал мне "Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up"

Код: Выделить всё

> c:=(m1sh,m1)->(m1sh/m1*cH2O+(m1-m1sh)/m1*cair):


> dmtk1:=(m1sh,m1,T1)->Thet(m1sh(t)+V1*muH2O*(alpha1*T1(t)^(betta1-1)-gamma1/T1(t))/R1)*(diff(m1sh(t),t)+V1*muH2O*(gamma1*diff(T1(t),t)/T1(t)^2+alpha1*T1(t)^(betta1-2)*diff(T1(t),t)*(betta1-1))/R1):

> f1:=diff(m1sh(t),t)=j0-j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))*m1sh(t)/m1(t)+j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))*m2sh(t)/m2(t)-dmtk1(m1sh,m1,T1):

> dmtk2:=(m2sh,m2,T2)->Thet(m2sh(t)+V2*muH2O*(alpha1*T2(t)^(betta1-1)-gamma1/T2(t))/R1)*(diff(m2sh(t),t)+V2*muH2O*(gamma1*diff(T2(t),t)/T2(t)^2+alpha1*T2(t)^(betta1-2)*diff(T2(t),t)*(betta1-1))/R1):

> f2:=diff(m2sh(t),t)=j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))*m1sh(t)/m1(t)-j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))*m2sh(t)/m2(t)-dmtk2(m2sh,m2,T2):
> f3:=diff(m1(t),t)=j0-j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))+j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))-dmtk1(m1sh,m1,T1):
> f4:=diff(m2(t),t)=j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))-j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))-dmtk2(m2sh,m2,T2):
> f5:=diff(c(m1sh(t),m1(t))*m1(t)*T1(t),t)=j0*cH2O*T0-c(m1sh(t),m1(t))*T1(t)*j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))+c(m2sh(t),m2(t))*T2(t)*j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))+dmtk1(m1sh,m1,T1)*(lambda1-cH2O*T1(t)):
> f6:=diff(c(m2sh(t),m2(t))*m2(t)*T2(t),t)=c(m1sh(t),m1(t))*T1(t)*j1(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))-c(m2sh(t),m2(t))*T2(t)*j2(m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t))+dmtk1(m2sh,m2,T2)*(lambda1-cH2O*T2(t)):
> muH2O:=0.018: muair:=0.029:R1:=8.3:V1:=1000:V2:=60000:gamma1:=10^(-20):alpha1:=6*10^(-
> 24):betta1:=11:cair:=1000:cH2O:=2000:A0:=1:T0:=400:S1:=2:j0:=3:g0:=9.8:lambda1:=2000000:
> Thet:=(x) -> Heaviside(x):F1:=(x) -> abs(x)*Heaviside(x):
> #Thet:=(x) -> `if`(x<0,0,1):F1:=(x) -> `if`(x<0,0,x):
> #Thet:=(x) -> 1/(exp(-50*x)+1):F1:=(x) -> abs(x)*Thet(x):
> deltaP2:=(m1sh,m2sh,m1,m2,T1,T2)->2*(m2/V2)*S1^2*R1*((m1sh/muH2O+(m1-m1sh)/muair)/V1*T1-
> (m2sh/muH2O+(m2-m2sh)/muair)/V2*T2):
> deltaP1:=(m1sh,m2sh,m1,m2,T1,T2)->2*(m1/V1)*S1^2*R1*((m1sh/muH2O+(m1-m1sh)/muair)/V1*T1-
> (m2sh/muH2O+(m2-m2sh)/muair)/V2*T2):
> j1:=(m1sh,m2sh,m1,m2,T1,T2)->sqrt(deltaP1(m1sh,m2sh,m1,m2,T1,T2)+2*m1/V1*A0*(m2/V2-m1/V1)
> *g0*S1^(5/2)):
> j2:=(m1sh,m2sh,m1,m2,T1,T2)->sqrt(F1(-deltaP2(m1sh,m2sh,m1,m2,T1,T2)+2*m1/V1*A0*(m2/V2-
> m1/V1)*g0*S1^(5/2))):
> j1(m1sh,m2sh,m1,m2,T1,T2):
> sys:=f1,f2,f3,f4,f5,f6:
> func:={m1sh(t),m2sh(t),m1(t),m2(t),T1(t),T2(t)}:
> m0:=m1sh(0)=0,m2sh(0)=0,m1(0)=1200,m2(0)=72000,T1(0)=350,T2(0)=350:
> ans:=dsolve({sys,m0},func,numeric,range=0..1);


https://www.dropbox.com/s/zugoypfg37dz5sz/may_be_no_trace.xmcdz
https://www.dropbox.com/s/coedmxbv5tyh75h/may_be.xmcdz
https://www.dropbox.com/s/iegrz5zk4ng5xme/third_try.mws

P.S. у меня вроде бы такая же версия Mathcad из того же источника. Версия: 15.0 (M010 [MC15_M010_20110622])

[Maxakch]

Я выше ошибся, не bp(), а pause(). bp() - так называлась моя аналогичная отладочная функция для моего отладчика.

Я нашел в интернете, упоминание об отладчике написанного uni, значит это вы. http://moryakin.tk/priboristu/otladchik_dlya_mathcad.html

П.С. Очень редко кто использует вывод отладочной информации таким образом. Такой метод использования отладочных возможностей Mathcad нужно брать на вооружение. Вряд ли даже на форуме PTC где-то можно увидеть такое.

да уж, я такого нигде не видел. Я да же не представлял, что в MathCad можно такие фокусы вытворять.
Такие вещи просто необходимо в FAQ прикреплять.

[uni]

После того как получены "живые" графики становится уже проще, т.к. наглядно виден процесс работы решателя. Далее уже всё зависит от теории и про "физичность" тут нужно думать самостоятельно.

По поводу Maple варианта нужно обратиться на соседний форум.

Я когда-то давно написал свой отладчик, но он был рассчитан на использования некоторых недокументированных приёмов, которые после MC11 убрали. Была возможность передавать в пользовательские библиотеки вложенные массивы с любыми типами данных.

Вот как работа моего отладчика выглядит на вашем примере. Там есть несколько окон: текущие значения переменных, журнал отладки и специальные окна для визуализации матриц. К сожалению, возможности отладчика не столь велики по сравнению со стандартными средствами, которые появились в MC13.



Вот, к примеру, что написано про отладку в современной книжке по MC15 "Инженерные расчёты в Mathcad 15. Учебный курс, 2011 г.":

Код: Выделить всё

Примечание (см. стр. 13): В Mathcad 13 появились средства отладки программ: панель Debug, функции trace и pause, - но реальной помощи в отладке программ они не оказывают.


Это всё на что хватило опыта и воображения у достаточно продвинутого автора, нашего коллеги, Евгения Макарова (книжка есть на трекере nnm-club.ru).

Причину тому я могу дать. Как и Поршнев, так и Макаров имеют смутное представление о программировании на ЯВУ. Это видно по их нелепым опечаткам в книгах. Я уже второй раз наблюдаю как в книжках по Mathcad пользовательские dll обзывают DDL-ками. Что показывает уровень владения материалом (см. стр. 29 "Панель контроля" - там вообще какая-то ерунда написана, что удивляет).

А если бы уважаемые авторы изучали Mathcad более полно, в том числе умели бы создавать эти самые dll, то им могло бы прийти в голову использовать некоторые приёмы из ЯВУ в самом Mathcad, как это сделал я. Но, увы, имеем книжку 2011 г., а материал по отладке там школьного уровня copy-paste (да-да, именно так автор рекомендует отлаживать программы - разборкой и сборкой, пока не заработает).

[Maxakch]

После того как получены "живые" графики становится уже проще, т.к. наглядно виден процесс работы решателя. Далее уже всё зависит от теории и про "физичность" тут нужно думать самостоятельно.

Большое спасибо за помощь. Буду думать, возможно, что то действительно необходимо поменять в теории. В ней был сделан ряд достаточно нахальных упрощений.

Да, я бы с удовольствием почитал Вашу книжку по Mathcad. Уверен, что у Вас набралась куча материалов по этому вычислительному пакету.

P.S. как и ожидалось, расчет с помощью rkfixed дал такой же результат как и Radau при увеличении количества шагов первого.