Чем нынче такое делать?

> Как это делают нормальные люди?

Строят линии уровня 0 прямо в гнуплоте.

octave

Во программа, кстати...

#include <iostream>
#include <cmath>
using namespace std;

const double t1=-5, t2=2.34, g=0.03, b=50, x_tick=1, z_tick=0.5, max_err=0.0001;
const int max_x=250, max_z=250;
const double t_needed=0;

double t(double x, double z)
{
  return( t1+(t2-t1)/M_PI*(atan((b+x)/z)+atan((b-x)/z))+g*z );
}

void say_x_z(double x, double z)
{
  cout << x << " " << z << endl;
}

double err(double x, double z)
{
  return(fabs(t(x,z)-t_needed));
}

int main()
{
  for(double z=0 ; z<=max_z ; z=z+z_tick)
  {
    if(z==0) {z=z_tick;}	//z не может быть равным нулю, т.к. функция при таком z не определена
    double t_last=t(0,z);
    for(double x=x_tick ; x<=max_x ; x=x+x_tick)
    {
      double t_curr=t(x,z);
      if((t_curr<=t_needed & t_last>=t_needed)|(t_curr>=t_needed & t_curr<=t_needed))
      {
        double delta_x=x_tick/2;
        x=x-delta_x;
        while(err(x,z)>max_err)
        {
          delta_x=delta_x/2;
          if(err(x+delta_x,z)>err(x-delta_x,z))
          {
            x=x-delta_x;
          }
          else
          {
            x=x+delta_x;
          }
        }
        say_x_z(x,z);
      }
      t_last=t_curr;
    }
  }
}

Тут вместо a(x,y) - t(x,z). Если лень читать... алгоритм таков: идём по Z с постоянным шагом. Для каждого Z идём по оси X до тех пор, пока t(x,z) не сменит знак. Когда t(x,z) меняет знак, методом половинного деления — точно не знаю, как этот (велосипедный?) метод называется — находим такой x между иксом до и после смены знака, чтобы ошибка не превышала допустимую.

не подкинешь толковый мануал?

> Строят линии уровня 0 прямо в гнуплоте.

Это как?

да, gnuplot и правда умеет строить изолинии o_O

чёрт.

Да, я это уже нагуглил (http://sovety.blogspot.com/2008/02/gnuplot-gri-pylab.html), но метод плох тем, что нельзя задать погрешность (можно задать только частоту семплирования)

>не подкинешь толковый мануал?

любой мануал по Matlab с небольшими изменениями. Довольно неплохая документация на официальном сайте.

man Graph Implicit Function

man Interval arithmetic

Задача является крайне сложной, так как решением может являться отдельная точка, кривая или область, а также их произвольное объединение.

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

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

Даже имеется готовые программы, которая все это реализуют. Я по крайней мере знаю вот эту:

http://www.peda.com/grafeq/

У меня банальная задача по геокриологии, тут по мерзлоте течет речка. Температура речки положительна, температура грунта на поверхности — отрицательна, плюс имеется теплопоток из недр. Очевидно, что если бы речки не было, то нулевая изотерма (линия, где t=0) была бы параллельна поверхности. Речка же тупо «выгибает» эту изотерму наверх , а в моём случае создаёт ещё и разрыв мерзлоты, в общем нулевая изотерма одна + есть небольшой диапазон глубин, где нуля нет вообще. Так что усложнять не надо.