Как работает Find в Mathcad

Здравствуйте.
Нужно написать програмку по вычислению некоторых формул.
Исходники и пример приводятся в маткаде.
Вроде бы проблем с большинством формул нету но попалась функция Find и я даже не могу толком сформулировать запрос чтобы решить вот такую хрень:

image

По отдельности все интегралы беруться отлично. Но в данном случае как я понял из мануалов маткада функция Find находит корень решения Lx1 при котором уравнение равно нулю.

Все бы ничего но я понятия не имею как вывести неизвестную из под интеграла… нас такому даже в универе не учили. Где взять исходник этой функции??
Может кто нибудь помочь??

Как вариант:
гугл подсказывает, что аналог в Matlab это fsolve —> он есть и в Octave —> исходники где-то тут (там в других файлах несколько реализаций для разных типов матриц), или через C++ API или stdin/stdout обертки

Блин … это дичь какая то … (((((

Тут функция в радикалах не берётся. Это значит что она решается только численно. А в Макаде тупо ошибки и он неверно считает.

1 лайк

Интеграл заменяется на функцию erf(x)
integral e^(-t^2/2) dt = sqrt(π/2) erf(t/sqrt(2)) + constant
Подставляются приделы интегрирования.
integral \limits ^{Lxt/2} _0 e^(-t^2/2) dt =sqrt(π/2) erf(Lxt/(2sqrt(2))) -0

Функция гладкая поэтому решается методом уточнения корня(Методом Ньютона).
erf(x) считается как интеграл по методу трапеций.
Или можно из алглиба выдрать приближённую формулу
https://www.alglib.net/specialfunctions/distributions/normal.php

Слабо верится в ошибки маткада. А вы ваш вариант можете посчитать??
А что за функция erf?? это откуда?

Есть и в маткаде, и в матлабе/октаве.

А в математике?

Это из теории вероятности и метрологии. Честно уже плохо помню что там к чему.

Возьмем штангенциркуль и нарисуем круг и начнём измерять его линейкой(на телефоне поставить программу рулетка). Фактический размер 6 мм. Но каждый раз у нас есть ошибки измерения 6,0 то 5,9 то 5,8. Математическое ожидание каждый раз изменяется от числа измерений. Ясно что размер круга постоянный, а ошибки связаны с методом измерения. Ошибки не могут иметь произвольные значения и имеют интервал допустимости. А так же распределение ошибок имеют некоторый вид erf(x) . Error function сокращенно erf.

Если плотность вероятности случайной величины задана функцией Гаусса, а мы хотим найти вероятность события, то мы интегрируем плотность и получаем функцию ошибок(error function).
Исторически сложилось, что функция erf связана именно с Гаусса, а не с любой другой функцией.
Вот так вот выглядит определеннее функции ошибок:
CodeCogsEqn (1)

Ньютона делать не стал сделала простым перебором. Поэтому считает около 10 секунд.

function erf(x:Real):Real;
var t,dt,h1,h2,s:Real;
 i:Integer;
begin
dt:=0.0001;
s:=0;
h1:=0;
for i:=0 to Round(x/dt)-1 do
  begin
  t:=i*dt;
  h2:=exp(-t*t);
  S:=S+(h1+h2)/2*dt;
  h1:=h2;
  end;
result:=2/sqrt(Pi)*s;
end;

procedure TForm1.Button1Click(Sender: TObject);
var Lx1,Lx,sigma_ix,delta_x,y:Real;
 i:Integer;
 y0,Lx0:Real;
begin
Lx:= 0.25;
sigma_ix:=0.125;
delta_x:= 0.048;
for i:=1 to 10000 do
  begin
  Lx1:=i/1000;
  y:=1/Lx1*erf(Lx1/2/sqrt(2))-1/Lx*erf(1/sqrt(2)*Lx/(2*sqrt(1+(sigma_ix+delta_x)*(sigma_ix+delta_x))));
  Chart1.Series[0].AddXY(Lx1,y);
  if i=1 then
     begin
     y0:=y;
     Lx0:=Lx1;
     end;
  if abs(y-0)<abs(y0-0) then
     begin
     y0:=y;
     Lx0:=Lx1;
     end;
  end;

  Memo1.Lines.Add('Результат find');
  Memo1.Lines.Add(Format('y0=%.03f;  Lx1=%.03f;',[y0,Lx0]));

end;

PS. Были ошибки в коде поправил.

2 лайка

Ох ничеси… Респект. :+1: :+1: :+1:

PS: Странно, Вроде сделал все также но у меня почему то результат другой: 0.6658

        float GetS(float Lx1, float lx, float sigma_ix, float delta_x)
        {
            int i;
            float y0 = 0, Lx0 = 0, y = 0;

            //Lx:= 0.25;
            //sigma_ix:= 0.125;
            //delta_x:= 0.048;
            int MxSteps = 10000;

            for (i = 1; i < MxSteps; i++)
            {
                Lx1 = (float)i / (float)MxSteps;
                y = 1.0f / Lx1 * erf(Lx1 / 2.0f / Math.Sqrt(2.0)) - 1.0f / lx * erf(1.0f / Math.Sqrt(2) * lx / (2.0f * Math.Sqrt(1.0f + (sigma_ix + delta_x) * (sigma_ix + delta_x))));

                if (i == 1)
                {
                    y0 = y;
                    Lx0 = Lx1;
                }

                double dtx = Math.Abs(y0 - y);

                if (Math.Abs(y - 0) < Math.Abs(y0 - 0))
                {
                    y0 = y;
                    Lx0 = Lx1;
                }
            }

            return Lx0;
        }

        float erf(double x)
        {
            double t, dt, h1, h2, S;
            int i;
            dt = 0.0001;
            S = 0;
            h1 = 0;
            for (i = 0; i < Math.Round(x / dt) - 1; i++)
            {
                t = i * dt;
                h2 = Math.Exp(-t * t);
                S = S + (h1 + h2) / 2.0 * dt;
                h1 = h2;
            }
            return (float)(2.0 / Math.Sqrt(Math.PI) * S);
        }

Тут минус 1 лишняя. Это моя ошибка.

Ваше

Тут ошибка, у меня там в 10 раз больше шаг. Я границы для Lx1 проставлял жестко от 0 до 10 с шагом 0,001. Что-бы график построить и посмотреть что там два корня. В маткаде условии начинают с 0.25 что-бы отсечь первый корень и взять второй. Соответственно в Вашем коде это надо в параметры функции выносить.

 Lx1 = (float)i / (float)MxSteps*10;

Мы не можем искать корень с той же точностью, что и считаем интеграл. Поэтому интеграл я считаю на 10 точнее. И плюс я считаю на double вы на float. Думаю дело в типах и там теряется точность.

Если интеграл заменить рядом Тейлора будет по быстрее и поточнее.
image

Удалось добится нужного результата при таких вот параметрах:

 double GetS(double Lx1, double lx, double sigma_ix, double delta_x)
        {
            int i;
            double y0 = 0, Lx0 = 0, y = 0;

            //Lx:= 0.25;
            //sigma_ix:= 0.125;
            //delta_x:= 0.048;
            **int MxSteps = 100000;**  <------ такая точность
            for (i = 1; i < MxSteps; i++)
            {
                Lx1 = (double)i / MxSteps;
                y = 1.0f / Lx1 * erf(Lx1 / 2.0f / Math.Sqrt(2.0)) - 1.0f / lx * erf(1.0f / Math.Sqrt(2.0) * lx / (2.0f * Math.Sqrt(1.0f + (sigma_ix + delta_x) * (sigma_ix + delta_x))));

                if (i == 1)
                {
                    y0 = y;
                    Lx0 = Lx1;
                }
                
                if (Math.Abs(y - 0.0) < Math.Abs(y0 - 0.0))
                {
                    y0 = y;
                    Lx0 = Lx1;
                }
            }

            return Lx0;
        }

        double erf(double x)
        {
            double t, dt, h1, h2, S;
            int i;
            **dt = 0.00001;**  <--------такая точность
            S = 0;
            h1 = 0;
            for (i = 0; i < Math.Round(x / dt); i++)
            {
                t = i * dt;
                h2 = Math.Exp(-t * t);
                S = S + (h1 + h2) / 2.0 * dt;
                h1 = h2;
            }
            return (double)(2.0 / Math.Sqrt(Math.PI) * S);
        }

Считает конечно долго. Тейлора попробую завтра.
Огромное спасибо. )))

2 сообщения было перенесено в новую тему: А почему он жирным не выделил запись (внутри кода)?

Здравствуйте еще раз. Что то с формулой где то не порядок. Не пойму почему.
Сделал реализацию через ряд тейлора. Действительно стало шустрее считать.

double erfTeilor(double x)
        {
            var p1 = 2.0 / Math.Sqrt(Math.PI);
            var p2 = erfS(x, 0, 0);
            return p1 * p2;
        }

        double erfS(double x, double prevX, int n)
        {
            double chisl = Math.Pow(-1.0, n) * Math.Pow(x, (2 * n) + 1);
            double znam = (double)Fact((UInt64)n) * ((2.0 * n) + 1);
            double S = (chisl / znam);

            //if (prevX != 0 && Math.Abs(prevX - x) < 0.0001)
            if (n > 30)
            {
                return S;
            }
            else return S + erfS(x, S, n + 1);
        }

При вот таких параметрах дает правильный результат который сходится с маткадом:


Далее нужно считать для Lz1 и тут какая то дичь идет.
Вот что в маткаде:

Вот что у меня получается:

Где то потерялось 0,101

Вот реализация вычислений по вашему коду.

        double GetS(double L_, double l_, double sigma_ix, double delta_x)
        {
            int i;
            double y0 = 0, Lx0 = L_, y = 0;
            //Lx:= 0.25;
            //sigma_ix:= 0.125;
            //delta_x:= 0.048;
            int MxSteps = 100000;
            for (i = 1; i < MxSteps; i++)
            {
                L_ = ((double)i / MxSteps);               

                y = 1.0 / L_ * erfTeilor(L_ / 2.0 / Math.Sqrt(2.0)) - 1.0 / l_ * erfTeilor(1.0 / Math.Sqrt(2.0) * l_ / (2.0 * Math.Sqrt(1.0 + (sigma_ix + delta_x) * (sigma_ix + delta_x))));

                if (i == 1)
                {
                    y0 = y;
                    Lx0 = L_;
                }

                if (Math.Abs(y - 0.0) < Math.Abs(y0 - 0.0))
                {
                    y0 = y;
                    Lx0 = L_;
                }
            }
            return Lx0;
        }

Я до конца так и не понял всех преобразований до функции erf(x), но может быть где то закралась ошибка?? Врядли эти погрешности высплывают из за типов чисел.
Уже перепроверил и в математике и в коде … даже не знаю что поправить.
А как выяснилось это ключевой момент всех вычислений. Если тут дать ошибку то дальше все плывет … :sleepy: :sleepy: :sleepy: :sleepy: :sleepy:

У вас факториниал пополняется ограничте n>15

Да не …
Проблема вот тут:
При таком подходе результата выше 1 не получить.

 int MxSteps = 100000;
            for (i = 1; i < MxSteps; i++)
            {
                L_ = ((double)i / MxSteps); 

переделал вот так:

double GetS(double L_, double l_, double sigma_ix, double delta_x)
        {
            double y0 = 0, Lx0 = L_, y = 0;
            int MxSteps = 1000;
            double step = 0.001;
            L_ = step;
            for (double i = step; i < MxSteps; i += step)
            {
                L_ = i;
                y = 1.0 / L_ * erfTeilor(L_ / 2.0 / Math.Sqrt(2.0)) - 1.0 / l_ * erfTeilor(1.0 / Math.Sqrt(2.0) * l_ / (2.0 * Math.Sqrt(1.0 + Math.Pow((sigma_ix + delta_x), 2))));

                if (y < 0)
                {
                    Lx0 = Math.Round(i, round);
                    break;
                }
            }
            return Lx0;
        }

По хорошему надо контролировать уменьшение значения, но я контролирую только смену знака.
Работает идеально. Результаты один в один с маткадом.