Математические зарисовки ⇐ Другие темы о русской словесности и культуре

[rusak]

52.

У Рамануджана есть вторая формула приближенного расчета длины эллипса. Она более точная. На ее основе создано совсем уж замечательное выражение:



Чтобы не мучиться и не набивать текст формулы, даю его на здоровьице:

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

L= Pi*(a+b)*(1+3*(a-b)^2/((a+b)^2*(10+sqrt(4-3*(a-b)^2/(a+b)^2))))*(22/(7*Pi)-(22/(7*Pi)-1)*exp(-0.4896141612727e-5*ln(b/a)^(12.472927032747/ln(b/a)^.526617815392+2.356417700185)))

Точность такая при b/a от 1 до миллиона, что получилось диво дивное и чудо чудное. Сама же формула Рамануджана приемлема только для b/a < 20
Расчет четырех параметров выполнялся по программе:

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

open #1,"ell-dat.txt","r"
open #2,"ell-2ramanud-t.txt","w"
dim x(300),y(300),f(300)
z=.001:t=22/7/pi:for i=1 to 72:input #1 x(i),y(i):next i
a0= -1:b0= 1:c0=-1:d0=1:s1=10^100:nn=100000000
for j=1 to nn:a=a0*(1+z*(ran()-.5)):b=b0*(1+z*(ran()-.5)):c=c0*(1+z*10*(ran()-.5)):d=d0*(1+z*(ran()-.5))
s=0:for i=1 to 72:x=x(i):y=y(i)
f(i)=t-(t-1)*exp(a*log(x)^(b*log(x)^c+d))
s=s+(y-f(i))^2:next i:if s<=s1 then:print f,a,b,c,d,s:s1=s:a0=a:b0=b:c0=c:d0=d:f0=f:r0=r
a1=a:b1=b:c1=c:d1=d:sk=s:fi:next j:print a1,b1,c1,d1,sk

Файл данных "ell-dat.txt" - это 72 точки. Первая колонка - отношение b/a, вторая колонка - отношение точного периметра эллипса к периметру, рассчитанному по второй формуле Рамануджана. Привожу эти точки;

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

 1.1  1.000000000
 1.2  1.000000000
 1.3  0.9999999997
 1.4  1.000000000
 1.5  1.000000000
 1.6  1.00000000
 1.7  0.9999999998
 1.8  0.9999999999
 1.9  1.000000000
 2  0.9999999999
 3  1.000000033
 4  1.000000249
 5  1.000000851
 6  1.000001964
 7  1.000003629
 8  1.000005821
 9  1.000008487
 10  1.000011558
 20  1.000052396
 30  1.000093261
 40  1.000127453
 50  1.000155380
 60  1.000178317
 70  1.000197394
 80  1.000213474
 90  1.000227195
 100  1.000239034
 200  1.000304457
 300  1.000332172
 400  1.000347560
 500  1.000357376
 600  1.000364190
 700  1.000369203
 800  1.000373052
 900  1.000376096
 1000  1.000378564
 2000  1.000390105
 3000  1.000394128
 4000  1.000396174
 5000  1.000397416
 6000  1.000398248
 7000  1.000398845
 8000  1.000399299
 9000  1.000399652
 10000  1.000399932
 200000  1.000402370
 300000  1.000402414
 400000  1.000402434
 500000  1.000402448
 600000  1.000402456
 700000  1.000402462
 800000  1.000402467
 900000  1.000402472
 1000000  1.000402474
 2000000  1.000402486
 3000000  1.000402491
 4000000  1.000402493
 5000000  1.000402494
 6000000  1.000402495
 7000000  1.000402496
 8000000  1.000402496
 9000000  1.000402497
 10000000  1.000402497
 20000000  1.000402498
 30000000  1.000402499
 40000000  1.000402498
 50000000  1.000402499
 60000000  1.000402499
 70000000  1.000402499
 80000000  1.000402499
 90000000  1.000402500
 100000000  1.000402500

Было выполнено 100 млн циклов, при которых методом Монте-Карло удалось вычислить очень точно параметры формулы Рамануджана-Александрова. Сумма квадратов отклонений (внимание!) составила величину 1.4165 , умноженной на десять в минус двенадцатой степени!
[/color]

[rusak]

Итак, прошло почти три недели. Они целиком были посвящены поиску наилучшей аппроксимации периметра эллипса.
Нужно отметить важную деталь: если отталкиваться от учебных пособий и от статьи "Эллипс" в Википедии, то обычно горизонтальная полуось a больше вертикальной полуоси b. Поэтому в формуле должно фигурировать отношение a/b , а не наоборот, как я делал это в прежних "Зарисовках".
Теперь дальше: последняя формула удивительно точна, но в ней нужно отношение b/a заменить на a/b . Это отношение находится в интервале от 1 до бесконечности. Еще более точная формула была получена 28 марта 2017 г. путем отказа как от числа Пи в формуле Рамануджана, так и от логарифма в поправочном коэффициенте. Вот так она выглядит:



Видно, что точность формулы Рамануджана резко начинает падать уже при отношении a/b>10. В то время, как эта точность в моей формуле уменьшается при прриближении к бесконечности. Это хорошо видно из второго рисунка:



Если же эллипс вытянут в вертикальном направлении ( то есть b>a), то в моей формуле эти параметры нужно просто поменять местами.

Важно понять следующее: приведенная выше приближенная формула периметра эллипса на сегодняшний момент является самой точной из множества других.

Текст фррмулы для Maple выглядит не столь устрашающе, как визуальная дробь на рисунке. Привожу листинг для копирования и проверки:

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

L := (22/7*(a+b))*(1+3*(a-b)^2/((a+b)^2*(10+sqrt(4-3*(a-b)^2/(a+b)^2))))/(1+(22/(7*Pi)-1)*exp(-((a/b)^.14597-1)^(3.40744/((a/b)^.389097-1)^.582524+1.12787)));

На этом "Зарисовки" по эллипсу завершены. Готовлю материал для доклада на конференции.

Как решалась эта задача. Программа на языке Yabasic по определению пяти оптимальных параметров знаменателя формулы Александрова:

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

open #1,"k.txt","r"
open #2,"k-otvet-17-112.txt","w"
dim x(200),y(200),f(200)
z=.0001:r=(22/7/pi-1)
for i=1 to 112
input #1 x(i),y
y(i)=y/10^6
w=w+1
print w,x(i),y(i)
next i
print r

b0= 1:c0= -1:d0=1:f0= 1:k0= 1 
s1=10^100:nn=800000000
for j=1 to nn
b=b0*(1+z*(ran()-.5))
c=c0*(1+z*(ran()-.5))
d=d0*(1+z*(ran()-.5))
f=f0*(1+z*(ran()-.5))
k=k0*(1+z*(ran()-.5))
s=0
for i=17 to 112
x=x(i)
f(i)=r*exp(-(x^f-1)^(b*(x^k-1)^c+d))
s=s+(y(i)-f(i))^2
next i
if s<=s1 then
print b,c,d,f,k,s
s1=s
b0=b:c0=c:d0=d:f0=f:k0=k
bk=b:ck=c:dk=d:fk=f:kk=k:sk=s
fi
next j
print #2,bk,ck,dk,fk,kk,sk

Файл данных - это невязка измененной формулы Рамануджана по сравнению с точными значениями:

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

            	        1.  402.499
                      1.2  402.499 
                      1.4  402.499 
                      1.6  402.499 
                      1.8  402.499 
                      2.0  402.499 
                      2.2  402.498 
                      2.4  402.494 
                      2.6  402.491 
                      2.8  402.481 
                      3.0  402.466 
                      3.2  402.444 
                      3.4  402.413 
                      3.6  402.372 
                      3.8  402.318 
                      4.0  402.250 
                       5.  401.648 
                       7.  398.869 
                       9.  394.010 
                      11.  387.528 
                      13.  379.950 
                      15.  371.712 
                      17.  363.129 
                      19.  354.430 
                      21.  345.764 
                      23.  337.235 
                      25.  328.909 
                      27.  320.827 
                      29.  313.012 
                      35.  291.254 
                      40.  275.010 
                      45.  260.342 
                      50.  247.081 
                      55.  235.064 
                      60.  224.141 
                      65.  214.179 
                      70.  205.063 
                      75.  196.694 
                      80.  188.985 
                      85.  181.864 
                      90.  175.265 
                      95.  169.135 
                      100.  163.427 
                      110.  153.108 
                      120.  144.036 
                      130.  135.996 
                      140.  128.823 
                      150.  122.381 
                      160.  116.563 
                      170.  111.284 
                      180.  106.470 
                      190.  102.062 
                      200.  98.012 
                      210.  94.276 
                      220.  90.819 
                      230.  87.611 
                      240.  84.625 
                      250.  81.840 
                      260.  79.234 
                      270.  76.793 
                      280.  74.499 
                      290.  72.341 
                      300.  70.306 
                      350.  61.654 
                      400.  54.921 
                      450.  49.526 
                      500.  45.108 
                      550.  41.420 
                      600.  38.294 
                      650.  35.610 
                      700.  33.281 
                      750.  31.240 
                      800.  29.437 
                      850.  27.832 
                      900.  26.394 
                      950.  25.099 
                      1000.  23.925 
                      1500.  16.316 
                      2000.  12.389 
                      2500.  9.989 
                      3000.  8.370 
                      3500.  7.204 
                      4000.  6.322 
                      4500.  5.634 
                      5000.  5.081 
                      5500.  4.627 
                      6000.  4.247 
                      6500.  3.926 
                      7000.  3.649 
                      7500.  3.409 
                      8000.  3.198 
                      8500.  3.012 
                      9000.  2.847 
                      9500.  2.699 
                      10000.  2.566 
                      15000.  1.716 
                      20000.  1.290 
                      25000.  1.033 
                      30000.  0.862 
                      35000.  0.740 
                      40000.  0.647 
                      45000.  0.576 
                      50000.  0.518 
                      55000.  0.471 
                      60000.  0.432 
                      65000.  0.398 
                      70000.  0.370 
                      75000.  0.346 
                      80000.  0.324 
                      85000.  0.305 
                      90000.  0.288 
                      95000.  0.273

Первая колонка - это отношение a/b при b=1. Вторая колонка - это (L0/R-1)*1000000 , где L0 - точное значение периметра эллипса; R - приближенное значение периметра по формуле Рамануджана, в которой число Пи заменено на его приближенное значение 22/7.
Рассчет длился чуть меньше трех суток.

[rusak]

Посоветовали умные люди оптимизировать структуру формулы приближенного расчета периметра эллипса и создать слайд для презентации. Уменьшенную копию формулы Александрова привожу здесь:

[rusak]

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

L := (22/7)*(a+b+3*(a-b)^2/(10*(a+b)+sqrt(a^2+14*a*b+b^2)))/(1+(22/(7*Pi)-1)*exp(-((a/b)^.14597-1)^(3.40744/((a/b)^.389097-1)^.582524+1.12787)));

[rusak]

Вчера и сегодня занимался анализом 17 формул расчета периметра эллипса. Не буду их тут рассматривать, а просто дам листинг для системы Maple. Кому интересно, всегда их выведет на свет божий:

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

restart;s := .825056176207: L1 := 4*(a+b)-(2*(4-Pi))*a*b/((1/2)*a^s+(1/2)*b^s)^(1/s);A := .518083: B := .538323: C := -.263801: D1 := 0.97480e-1: K := 1.001249: L2 := 4*a*((A+B*(b/a)^C)^(-1/D1)+K);L3 := Pi*sqrt(2*(a^2+b^2)-(1/2)*(a-b)^2);L4:=Pi*(3*(a+b)-sqrt((3*a+b)*(a+3*b)));L5 := Pi*(a+b)*(1+3*((a-b)/(a+b))^2/(10+sqrt(4-3*((a-b)/(a+b))^2)));x := log(2)/log((1/2)*Pi): L6 := 4*(a^x+b^x)^(1/x);L7 := Pi*sqrt(2*(a^2+b^2));L8 := Pi*(a+b);L9 := Pi*(3/2*(a+b)-sqrt(a*b));L10 := 2*Pi*a*sqrt(3/4+(1/4)*b/a-3/8*(1-(b/a)^2));L11 := a*(2.28*(b/a)^1.44+4);L12 := b*(2.5885*(a/b)^(1.2417*(b/a)^0.3155e-1)+3.6947);L13 := Pi*sqrt(2*(a^2+b^2)-(1/4)*(a-b)^2);L14 := (1/2)*Pi*a*(3*(1+b/a)-2*sqrt(b/a));L15 := Pi*(a+b)*(1+3*(a-b)^2/((a+b)^2*(10+sqrt(4-3*(a-b)^2/(a+b)^2))))*(1+(22/(7*Pi)-1)*((a-b)/a)^33.697);L16 := (22/7)*(a+b+3*(a-b)^2/(10*a+10*b+sqrt(a^2+14*a*b+b^2)))/(1+(22/(7*Pi)-1)*exp(-((a/b)^.14597-1)^(3.40744/((a/b)^.389097-1)^.582524+1.12787)));L17 := 4*b*EllipticE(sqrt(-a^2/b^2+1));

Формулы пронумерованы от L1 до L17. Отмечу только шесть формул-рекордсменов по точности вычислений. Это чтобы была понятна последующая таблица. Итак:
L1 - формула Кантрелла
L4 - первая формула Рамануджана
L5 - вторая формула Рамануджана
L15 - формула из японской Википедии (статья "Эллипс")
L16 - формула Александрова
L17 - точная формула через эллиптический интеграл.
Теперь сопоставительная таблица:


Красная черта - первое место по точности среди приближенных формул, зеленая черта - 2-е место, синяя черта - 3-е место.
Видно, что в основном побеждает моя формула. Формула японцев лишь локально мощно вырвала первенство в точности вычислений. Меня это удивило. Решил проверить это на графике и понять. Вот как меняется точность вычисления функций L15 и L16 по сравнению с абсолютом L17:



Видно, что японская формула L15 пересекает ось Оа как раз при а=10. Так уж случайно получилось. Но хорошо видно, что моя аппроксимация имеет значительно меньшую амплитуду отклонения от точного решения.

Как и ожидалось, формулы Рамануджана сдают свои позиции по мере увеличения разницы между значениями полуосей эллипса. Зато при очень большой разнице на второе место выходит формула Кантрелла.
David Cantrell - современный математик (родился в 1949 г.), прекрасно владеющий теорией эллиптических функций, программист, специалист по аппроксимациям. Это вам не Рамануджан, который всего добивался чудовищной интуицией, карандашом и блокнотом. Кантрелл пользовался всем арсеналом электронных математических систем и всего лишь лет 15 назад вывел свою формулу (если хорошо погуглить, то можно найти его статью - начинать с http://www.numericana.com/fame/index.htm и http://www.numericana.com/answer/ellipse.htm#cantrell ). Вот с таким динозавром пришлось мне соперничать. Сопоставительный график наглядно показывает, на сколько порядков точность по модели Кантрелла меньше моей ( синяя кривая Александрова буквально ползет по оси абсцисс):



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

[rusak]

Времени нет, завтра лететь в Дубай. Прорыв такой в части периметра эллипса, что не хватит и ста страниц. Довел до совершенства уже две формулы.

Финишный итог:
ПЕРИМЕТР ЭЛЛИПСА
первая формула Александрова
вторая формула Александрова:



Первая формула (менее точная - от первой формулы Рамануджана):

L:= 4*(3*(a+b)-sqrt((3*a+b)*(a+3*b)))/((3-sqrt(3))*(1+(4/(Pi*(3-sqrt(3)))-1)*exp(-251*((a/b)^0.00869-1)^(1.23/((a/b)^1.15-1)^.527+1.48))));

Вторая формула (от второй формулы Рамануджана):

L:= (22/7)*(a+b)*(1+3*(a-b)^2/((a+b)^2*(10+sqrt(4-3*(a-b)^2/(a+b)^2))))/(1+(22/(7*Pi)-1)*exp(-157*((a/b)^0.012-1)^(2.69/((a/b)^.929-1)^.588+1.58)));

[rusak]

Целый месяц убил на обновление рекорда точности расчета периметра эллипса. Одной формулой сделать это не удалось (то есть формула Александрова L2 так и остается лучшей в мире), но двумя аппроксимациями можно на порядок повысить эту самую капризную точность. Помог мне коллега из форума dxdy.ru под ником Xaositect. Он сделал потрясающе точное разложение для значений a/b > 13. Для интервала же от 1 до 13 нужно было дать иную аппроксимацию. И тут пригодилась моя формула, но с иными оптимизированными шестью параметрами! До значения 13 точность безукоризненная. Рамануджан уже отдыхает окончательно. Все, что я сказал, наглядно демонстрирует коллаж:



Все познается в сравнении. Поэтому сравним мою лучшую формулу L2 с формулой, которая оптимально аппроксимирует периметр эллипса для интервала значений a/b от 1 до 13:



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

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

open #1,"k.txt","r"
open #2,"21с.txt","w"
dim x(200),y(200),f(200)
z1=10:z=10:r=(22/7/pi-1)
zz=.5
for i=1 to 176
input #1 x(i),y
y(i)=y/10^6
w=w+1
print w,x(i),y(i)
next i
nn=1300000
for v=1 to 500
s1=10^100
print v
w0=120:c0=77:d0=327:f0=238:k0=275:t0=141
for j=1 to nn
w=int(w0+z1*ran()-z*.3)
c=int(c0+z1*ran()-z*zz)
d=int(d0+z1*ran()-z*zz)
f=int(f0+z1*ran()-z*zz)
k=int(k0+z1*ran()-z*zz)
t=int(t0+z1*ran()-z*zz)
s=0
for i=2 to 21
x=x(i)
f(i)=r*exp(-w*(x^(c/10000)-1)^((d/100)*(x^(f/100)-1)^(-(k/1000))+(t/100)))
s=s+(y(i)-f(i))^2
next i
if s<=s1 then
print w,c,d,f,k,t,s
s1=s
w0=w:c0=c:d0=d:f0=f:k0=k:t0=t
wk=w:ck=c:dk=d:fk=f:kk=k:tk=t:sk=s
fi
next j
print wk,ck,dk,fk,kk,tk,sk
print #2,wk,ck,dk,fk,kk,tk,sk
next v

Это моя одна из самых уникальных программ. Дело в том, что если комбинаторно прокручивать все шесть параметров формулы, то потребуются годы непрерывного счета. Здесь же я придумал способ случайного ветвления. Представим себе дерево, по стволу которого ползет муравей. Случайным образом он сворачивает на одну из веток, далее случайно поворачивает на любой отросток и так далее, пока не окажется на самом конце ветки. Этот конец ветки - есть некий локальный оптимум. Чтобы вручную не запускать прогу, я предусмотрел автоматическую реализацию пятисот совершенно независимых прогонов. И затем в компактной распечатке ищу глобальный оптимум. Одна ночь - это как раз 500 локальных минимумов. Несколько ночей потребовалось, чтобы выявить глобальный минимум, и эти коэффициенты видны на первом рисунке данного поста.

Теперь можно и отдохнуть.

[rusak]

Периметр эллипса продолжаю улучшать. Аппроксимирующая программа крутится уже более месяца - настолько сложной является кривая.
Но два последних дня посвятил удивительно красивой школьной задаче. Оказалось, что не всё еще известно о фигуре трапеция. На математическом форуме был задан вопрос: если известны боковые стороны и диагонали трапеции, то как можно найти ее основания? Потуги мои привели к новому шедевру. Вот окончательная картинка:



Посмотрите на эти формулы - они потрясающе гармоничные и никакой фантазии не хватит, чтобы их изобрести. Это - поэзия алгебры!
В системе Maple написал программу проверки правильности вычислений диагоналей и оснований трапеции, а также проверку существования фигуры. Если в конце расчетов появится сообщение "no!", то это будет означать, что нарушено хотя бы одно из четырех ограничений. Итак, текст проги:

restart; a0 := 3.; b0 := 6.; c := 5.; d := 4.; d1 := sqrt(a0*b0+d^2+b0*(c^2-d^2)/(b0-a0)); d2 := sqrt(a0*b0+c^2-b0*(c^2-d^2)/(b0-a0)); a := sqrt(((c^2-d1^2)^2-(d^2-d2^2)^2)/(2*(c^2-d^2+d1^2-d2^2))): b := sqrt(((c^2-d2^2)^2-(d^2-d1^2)^2)/(2*(c^2-d^2-d1^2+d2^2))): a := evalf(a); b := evalf(b); print(d1, a+c, b+d); print(d2, a+d, c+b); if d2 >= a+d or d2 >= c+b or d1 >= a+c or d1 >= b+d then print("no!") end if

Жду с нетерпением отзывы специалистов со всего белого света.

[rusak]

Опубликовал коротенькую статью в своем блоге
http://renuar911.blog.ru/?year=2017&month=07&day=8

Однако! Эти формулы удалось обнаружить во французской Википедии! В несколько ином, правда, виде. Но тождественном. Значит, у меня все верно!
https://fr.wikipedia.org/wiki/Trap%C3%A8ze

Недолго думая, поместил свои формулы в нашу, российскую часть Википедии.
В Maple составил прогу проверки двух формул. А также проверку на возможность самой трапеции. Если в конце распечатки появится возглас no!, значит, трапеция при исходных данных невозможна (будет наблюдаться нарушение одного из четырех ограничений: d1<a+c;d1<b+d; d2<a+d; d2<c+b ):

restart; a0 := 3.; b0 := 6.; c := 5.; d := 4.; d1 := sqrt(a0*b0+d^2+b0*(c^2-d^2)/(b0-a0)); d2 := sqrt(a0*b0+c^2-b0*(c^2-d^2)/(b0-a0)); a := sqrt(((c^2-d1^2)^2-(d^2-d2^2)^2)/(2*(c^2-d^2+d1^2-d2^2))): b := sqrt(((c^2-d2^2)^2-(d^2-d1^2)^2)/(2*(c^2-d^2-d1^2+d2^2))): a := evalf(a); b := evalf(b); print(d1, a+c, b+d); print(d2, a+d, c+b); if d2 >= a+d or d2 >= c+b or d1 >= a+c or d1 >= b+d then print("no!") end if

Произвел интересное кинематическое исследование шарнирной рамы, состоящей из элементов боковых сторон трапеции и диагоналей. Расписал общие формулы изменения высот точек B и C и нашел точку, которая определяет трапецию. Все расчеты и графики в системе Maple привожу тут же:


restart; yB := sqrt(1/2*(c^2+d2^2-1/2*(x^2+(c^2-d2^2)^2/x^2))); yC := sqrt(1/2*(d1^2+d^2-1/2*(x^2+(-d^2+d1^2)^2/x^2))); c := 50; d := 60; d1 := 60; d2 := 90; plot([yB, yC], x = -30 .. 130); x := 67.91; evalf(yB); evalf(yC)

Были опасения, что решений не одно, а несколько. После приведенного анализа сомнения развеялись и в данной зарисовке только что доказано: решение единственное.

[rusak]

Во время отдыха в Бенидорме ко мне на математическом форуме обратился студент с просьбой вывести формулу объема усеченного клина. И еще показать, как эта формула может трансформироваться в объем пирамиды. На пляже подобные задачки решаются легко, как по маслу. Результаты оформил уже поздно ночью:



Здесь на клин смотрим сверху. Основание - прямоугольник.
Это мои пять копеек в компактные шпаргалки для любителей геометрии.

[rusak]

Яндекс уже читает мой основной рисунок по трапеции! Он в левом верхнем углу - я его обвел красной линией:

[rusak]

Гугл тоже заметил мой рисунок! Хотя и постеснялся дать его покрупней:

[rusak]

Сегодня весь день посвятил биссектрисе. Аналитической биссектрисе. Циркулем и линейкой ее любой школьник построит. А вот уравнение написать - это слабо и сапожнику, и дояру. Короче, даю два рисунка. Вроде, все понятно и доступно даже поэтам.


Также и тут в неявном виде http://www.pm298.ru/reshenie/biss.php
А вот так мучаются все студенты: http://dokak.ru/matematika/135297-kak-n ... trisy.html



Поместил свое решение в Вики. Надеюсь, что не удалят.
Вчера отладил программу в системе Maple. Привожу текст:
with(plots): a1 := -1.5; b1 := 6; a2 := 4.5; b2 := 3; x0 := (b2-b1)/(a1-a2); y0 := a1*x0+b1; q1 := plot({a1*x+b1, a2*x+b2}, x = x0-1 .. x0+1, y = y0-1 .. y0+1, color = black, thickness = 3,scaling=CONSTRAINED): q2 := plot({(a1*sqrt(a2^2+1)-a2*sqrt(a1^2+1))*x/(sqrt(a2^2+1)-sqrt(a1^2+1))+(b1*sqrt(a2^2+1)-b2*sqrt(a1^2+1))/(sqrt(a2^2+1)-sqrt(a1^2+1)), (a1*sqrt(a2^2+1)+a2*sqrt(a1^2+1))*x/(sqrt(a2^2+1)+sqrt(a1^2+1))+(b1*sqrt(a2^2+1)+b2*sqrt(a1^2+1))/(sqrt(a2^2+1)+sqrt(a1^2+1))}, x = x0-1 .. x0+1, y = y0-1 .. y0+1, color = red,scaling=CONSTRAINED): display(q1, q2);
В данном случае график такой (красные линии - биссектрисы):

Маленький шедевр! Завтра в Китае расскажу и покажу студентам.

[rusak]

Самые лучшие математические открытия мне почему-то удается делать вдали от России. Наверное, потому что меньше думается о политике, а больше об окружностях, точках и углах. Вот и в Гонконге все мысли были брошены на решение простой геометрической, но очень сложной в аналитическом плане задачи. Имеем угол A, который образуется пересечением двух прямых. Уравнения прямых: y1=a1*x+b1 и y2=a2*x+b2. Внутри этого угла имеется точка P, которая принадлежит вписанной в угол окружности. Нужно найти координату центра окружности и ее радиус.
При помощи циркуля и линейки задание выполняется без проблем. Строим биссектрису угла, строим в произвольном месте вспомогательную вписанную окружность, строим луч AP и далее путем построения параллельных линий находим не одну, а две вписанные окружности, удовлетворяющие условиям задачи. Но это геометрически.
С аналитикой дело усложняется стократно. В литературе ничего путного найти не смог, за исключением частных числовых примеров. Эти примеры и решаются численными методами.
Но мне желательно получить как раз общий метод в виде алгебраических формул. Два дня пришлось плотно поработать в системе Maple, чтобы вывести формулы, их существенно упростить и составить программу визуализации получаемых результатов при задании четырех параметров прямых и двух координат точки P. Математическая модель и пример расчета уместились на одном листе. Его я и предлагаю на рассмотрение:

Программа в Maple:
restart; with(plots): k := 2; xp := 1.58*2; yp := .32*2; a1 := 389/172; b1 := 0; a2 := 0; b2 := 0; if k = 1 then ab := evalf((a1*sqrt(a2^2+1)-a2*sqrt(a1^2+1))/(sqrt(a2^2+1)-sqrt(a1^2+1))); bb := evalf((b1*sqrt(a2^2+1)-b2*sqrt(a1^2+1))/(sqrt(a2^2+1)-sqrt(a1^2+1))) end if; k := k; if k = 2 then ab := evalf((a1*sqrt(a2^2+1)+a2*sqrt(a1^2+1))/(sqrt(a2^2+1)+sqrt(a1^2+1))); bb := evalf((b1*sqrt(a2^2+1)+b2*sqrt(a1^2+1))/(sqrt(a2^2+1)+sqrt(a1^2+1))) end if; q := -(a1*ab+1)^2/(a1^2+1); t := 2*((a1-ab)*(b1-bb)/(a1^2+1)+ab*(yp-bb)+xp); v := (b1-bb)^2/(a1^2+1)-xp^2-(bb-yp)^2; x1 := evalf((-t+sqrt(-4*q*v+t^2))/(2*q)); x2 := evalf((-t-sqrt(-4*q*v+t^2))/(2*q)); yb1 := ab*x1+bb; yb2 := ab*x2+bb; r1 := sqrt((xp-x1)^2+(yp-yb1)^2); r2 := sqrt((xp-x2)^2+(yp-yb2)^2); q1 := plot({-sqrt(r1^2-(x-x1)^2)+yb1, sqrt(r1^2-(x-x1)^2)+yb1, -sqrt(r2^2-(x-x2)^2)+yb2, sqrt(r2^2-(x-x2)^2)+yb2, ab*x+bb}, x = -.5 .. 25, y = -15 .. 15, scaling = CONSTRAINED, color = black):q2:=plot({a1*x+b1, a2*x+b2}, x = -.5 .. 25, y = -15 .. 15, scaling = CONSTRAINED,thickness=3): data := [[xp, yp], [x1, yb1], [x2, yb2]]: q3 := pointplot(data, symbol = circle,symbolsize = 15, color = "Red"): display(q2, q1, q3);
В результате получим рисунок:

Уже в четверг, 27 июля 2017 г. по данной проблеме сделал доклад в Университете. Скучно студентам не было.

[rusak]

Гугл заметно продвинул мои формулы для расчета оснований трапеции по известным боковым сторонам и диагоналям. Рисунок уже на второй строке! Ключевае предложение: "Формула Александрова для расчета трапеции". Ссылку дает как раз на rusforus.ru, за что выражаю признательность нашему форуму. Значит, он более популярный в мире, нежели все математические русскоязычные форумы. Показываю скриншот:

[rusak]

Со мной связался по эл. почте профессор одного из университетов и попросил доказать верность формул для биссектрис, основываясь на Примере 17 по ссылке http://www.diary.ru/~eek/p56326549.htm
Вот этот пример:

Я привык видеть уравнения в явном виде. Поэтому условия задачи такие: заданы уравнения сторон угла:
y=3x/4+7/4 и y=-15x/12+1/12
уравнения биссектрис получились следующими:
y=x/8+6/7 и y=-8x-43/4
Мои формулы для биссектрис можно записать более компактно:


Программа в системе Maple:
with(plots); a1 := 3/4; b1 := 7/4; a2 := -5/12; b2 := 1/12; v := sqrt((a1^2+1)/(a2^2+1)); y1 := (a2*v+a1)*x/(1+v)+(b2*v+b1)/(1+v); y2 := (-a2*v+a1)*x/(1-v)+(-b2*v+b1)/(1-v); q1 := plot({a1*x+b1, a2*x+b2}, x = -9 .. 6, y = -4 .. 5, scaling = CONSTRAINED, color = black, thickness = 3): q2 := plot({(-a2*v+a1)*x/(1-v)+(-b2*v+b1)/(1-v), (a2*v+a1)*x/(1+v)+(b2*v+b1)/(1+v)}, x = -9 .. 6, y = -4 .. 5, scaling = CONSTRAINED, color = purple, thickness = 1): display(q1, q2)

По этой программе получим ответы Примера 17:


Как видим, все совпало, но мой подход оказался намного проще.
Профессор оценил метод и будет его популяризировать. Я конечно обеими руками "за".

[rusak]

Опять про трапецию. Не дает покоя эта прекрасная фигура. Вывел новые полезные формулы, когда известны основания a и b , а также обе диагонали. Выражения оказались настолько красивыми, что невозможно их здесь не привести:


Если все подытожить, то получим великолепнейший справочный лист:

[rusak]

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



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



Черные числа - это размеры, которые вычисляются по найденным мною формулам. Формулы в программах просматриваются и поэтому не буду их специально выставлять.
Первый вариант рассчитывается по программе:

restart; yB := sqrt(1/2*(c^2+d2^2-1/2*(b^2+(c^2-d2^2)^2/b^2))); yC := sqrt(1/2*(d^2+d1^2-1/2*(b^2+(d^2-d1^2)^2/b^2)));xB:=sqrt(c^2-yB^2); xC := sqrt(d1^2-yC^2); yA := 0; xA := 0; yD := 0; xD := b; c := 50; d := 60; d1 := 90; d2 := 70;b := sqrt((d*(d2^2-c^2)+c*(d1^2-d^2))/(c+d)): yB := yB; evalf(yB); xB :=xB; evalf(xB); yC := yC; evalf(yC); xC := xC; evalf(xC); xD := b; evalf(xD); a :=evalf(sqrt((xB-xC)^2+(yB-yC)^2));b:=evalf(b); A := [xA, yA]: B := [xB, yB]: C := [xC, yC]: D1 := [xD, yD]: AB := [A, B]: BC := [B, C]: CD := [C, D1]: AD := [A, D1]: BD := [B, D1]: AC := [A, C]: plot({AB, AC, AD, BC, BD, CD}, x = min(0, xB) .. max(b,xC), y = 0 .. max(yB, yC), color = black, scaling = CONSTRAINED, thickness = 3);

Второй вариант:

restart; yB := sqrt(1/2*(c^2+d2^2-1/2*(b^2+(c^2-d2^2)^2/b^2))); yC := sqrt(1/2*(d^2+d1^2-1/2*(b^2+(d^2-d1^2)^2/b^2)));xB:=k*sqrt(c^2-yB^2); xC := sqrt(d1^2-yC^2); a := sqrt(((d2^2-d^2)*c+d*(-c^2+d1^2))/(c+d));b := sqrt((c*(d1^2-d^2)+d*(d2^2-c^2))/(c+d));yA := 0; xA := 0; yD := 0; xD := b: c := 50; d := 60; d1 := 80; d2 := 70;if d1>=d2 then k:=1 else k:=-1 end if: yB :=evalf(yB); xB :=evalf(xB); yC :=evalf(yC); xC :=evalf(xC); xD :=evalf(b): a :=evalf(a);b:=evalf(b); A := [xA, yA]: B := [xB, yB]: C := [xC, yC]: D1 := [xD, yD]: AB := [A, B]: BC := [B, C]: CD := [C, D1]: AD := [A, D1]: BD := [B, D1]: AC := [A, C]: plot({AB, AC, AD, BC, BD, CD}, x = min(0, xB) .. max(b,xC,xB), y = 0 .. max(yB, yC), color = black, scaling = CONSTRAINED, thickness = 3);

[rusak]

Только что полностью завершил поставленную задачу: зная боковые стороны трапеции и ее диагонали, найти аналитически основания трапеции, среднюю линию и медиану. Если средняя линия m есть полусумма оснований, то формулу для медианы трапеции (обозначаю ее как n ) мне не удалось встретить ни в Википедии, ни в литературных источниках, ни на картинках и даже на лекциях по геометрии в ютубе. А формула эта оказалась настоящим шедевром! Ее упростить мне помогли коллеги на математическом форуме. Но ближе к делу.
Формулы и пример расчета:

[rusak]

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



где S - площадь трапеции.

Что тут важно отметить? Раньше я решал такие задачи: известны 4 параметра трапеции и по ним находились два неизвестных. Здесь же по известным только основаниям a и b удалось выудить ценную информацию о соотношениях площадей простейших элементов - треугольников. Новое знание позволит сильно упрощать решение многих геометрических проблем.

Вместе с тем, продолжил и менее ценные задачи по трапеции. Они касаются вывода формулы для длины медианы по известным четырем параметром. Одна из них, формула (5), наверняка получена мной впервые, ибо ее вывод оказался чрезвычайно сложным. Но красота покорила не только меня. Показываю вывод трех формул, который я неделю назад продемонстрировал на конференции:

Что такое трапеция? Она по сути - усеченный треугольник. В треугольнике одной из важнейших линий считается медиана EG. Очевидно, отрезок FG является оставшейся частью медианы треугольника AED. То есть его, этот отрезок, можно назвать медианой трапеции ABCD.
Но почему-то везде и всюду любят говорить только о средней линии трапеции m и совершенно не интересуются медианой. Хотя наверняка в литературе формула для FG получена, но в учебниках не распространилась.




Попытаемся вывести длину медианы трапеции.
Из рисунка хорошо видно, что сторона трапеции c в линейной пропорции увеличивается с уменьшением основания a путем его перемещения вверх. И наступит такой момент, когда a=0. Это произойдет в точке E. Составим пропорцию:



Сделал проверку для "широких" и "высоких" трапеций:

[rusak]

Еще в десятом классе я сильно увлекся трапецией и выводил различные формулы, которые не фигурировали в учебниках. Тогда я вел дневник и, к счастью, он сохранил мои творческие мысли и результаты. Медиану обозначил через n . После рисунка и формул идет примечательное: "За эти формулы наша Л.Е. поставила круглую пятерку". Инициалы Л.Е. - это наша прекрасная учительница по математике Карманова Людмила Ефремовна (17 лет назад я ей посвятил небольшую статью http://renuar911.narod.ru/chislo_e.htm ).

[rusak]

О, Боже! Трапеция никак меня не отпускает и удивляет новыми сказочными открытиями. Решил посмотреть, какая формула получится для высоты h , если вывести ее в зависимости от параметров оснований a и b, а также двух диагоналей. Начал как обычно вручную, но запутался и призвал на помощь Maple18 и Вольфрам Альфа. Одно из упрощений оказалось невероятно красивым, - комбинаторным! Запомнить формулу теперь можно сразу и навсегда. Вот, удивляйтесь вместе со мной:



Аналогичным образом высоту трапеции можно найти даже по трем известным параметрам трапеции:



Это вытекает из того, что дан жесткий треугольник со сторонами b, d, d1. Подобную формулу можно составить и для второго треугольника со сторонами b, c, d2 . Вместо d ставим d2 и вместо d1 ставим c:



Проверка двух формул ( данные взяты с произвольно начерченного графика):



И еще два аналогичных тождества для h можно вывести по известным сторонам треугольника a, c, d1 и треугольника a, d, d2. Это по силам сделать теперь любому старшекласснику.

Столько прекрасных общих формул заложено в трапеции! И все они запрятаны в урезанной числовой форме - в виде тысяч частных задач. Думаю, настало время систематизировать общие формулы, ибо они облегчат жизнь школьников, студентов, инженеров, аспирантов...

[rusak]

К такому утешительному выводу пришел вчера:

[Кеша]

К такому утешительному выводу пришел вчера:

Дело/100, Ас-Салям/100 гроб/100подошёл/100?ждёт/100ДЬЯВОЛ/100

[rusak]

Если бы только 100. Чтобы полностью аналитически разрешить общую трапецию, в которой выявлено 8 параметров ( a b c d d1 d2 h n), нужно создать 1680 схем, подобных такой:



Работы примерно на год. Книгу буду делать в электронном представлении. Причем формулы дам в виде, как набивают в системе Maple или в любом калькуляторе. Можно быстро скопировать и сразу численно вычислять параметры. Современно и безошибочно! Каждая формула - 1 строка. Итого в каждой схеме будет 5 строк. Пусть лист (электронный) содержит 40 строк. Тогда на каждой странице поместится 8 вариантов схем. Потребуется 1680/8=210 страниц. Многовато, конечно, но реально. Зато это будет первое в мире полное аналитическое решение трапеции.