Глава 2. В которой  внук Андрюша

 узнает, что такое аппроксимация

 

     В  “Интернете”  я часто посещаю студенческий форум и в разделе «Математика»  помогаю многим в решении различных задач по высшей математике.  Две недели назад  некий аспирант по имени  modul  обратился ко всем с просьбой подобрать функцию, которая хорошо бы описывала экспериментальные данные. Вот такую таблицу он приложил:

Несколько  модераторов  и студентов  предпринимали  попытки  найти хотя бы приблизительную зависимость,  но так и не смогли этого  сделать  без явных нарушений. Например,  у одной дамы  кривая не проходила через начало координат.  Я уже забыл про задачу, но сегодня Андрюша прибежал ко мне после завтрака и спросил:

- Деда, что такое  “апраксимация”?

- А откуда ты такое слово услышал?

- По телевизору дядя сказал.

- Нужно говорить  аппроксимация.  В этом слове два “п”  и после буквы “p” идет буква “о”.

- Да-да, “оп-проксимация”.

- Первая буква “а”  , а не “о”.

-Ну хватит, дедуля! Лучше скажи, что это такое.

- Это, внучек, одно из основных понятий математики.  В общем случае это замена одних математических объектов другими, более простыми и удобными.  Например, мы с тобой проделаем  какой-нибудь опыт и получим много-много точек наблюдений, которые зафиксируем на бумаге. И мы увидим, что эти точки  не кое-как разбросаны по координатной сетке, а выстраиваются в какую-то гладкую линию. Конечно же,  зададимся вопросом: а можно ли подобрать такую функцию,  чтобы опытные точки явились как бы частью кривой. Такой вот подбор функции и носит название  “аппроксимация”.

- Давай  попробуем сделать аппроксимацию!

- Андрюша! Это же очень сложно. Хотя  я  в бытность  собаку съел на этом.

- Какую собаку, дедушка?

- Это так обычно говорят, когда много занимаются одним и тем же делом.

- Раз ты собаку съел и много занимался, то расскажи! Я тоже хочу  много заниматься.

- Хорошо. Я вот вспомнил про одну задачу, в ней как раз были опытные точки. Сейчас я найду  тему, поднятую аспирантом  с именем modul … А, вот он, этот аспирант! Вот и таблица… Видишь, даны ни много, ни мало, 32 точки. И до сих пор несчастному будущему  ученому  никто толком так и не помог. Придется  нам с тобой  сделать  для него доброе дело.

- Давай, дедуля! Сделаем  доброе  дело.

- Прежде всего строим точечный график, чтобы убедиться, что с данным набором точек можно успешно работать. Я тебе буду диктовать, ты набивай команды.

- Поехали!

 

with(plots):X:=[.2,.4,.6,.8,1,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6,2.8,3,3.2,3.4,3.6,3.8,4,4.2,4.4,4.6,4.8,5,5.2,5.4,5.6,5.8,6,6.2,6.4]:Y:=[.02,.07,.17,.36,.55,.75,.96,1.11,1.21,1.29,1.28,1.24,1.18,1.07,.95,.85,.72,.59,.5,.4,.31,.26,.19,.14,.12,.08,.05,.05,.03,.02,.02,.01]:g:=plot([[X[i],Y[i]]$i=1..32],x=0..6.5,style=POINT,symbol=CIRCLE,color=blue):display(g);

 

 

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

 

 

Тут три независимых параметра.  Если графически рассмотреть сомножители, то получаются  примерно  такие две кривые:

 

 

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

- Кто такой Монте-Карло?

- Это город такой в Европе. Там много игровых автоматов. Сейчас, Андрюшенька, я попытаюсь вспомнить  программирование на языке Yabasic.  Следи внимательно за тем, как  набиваю текст  программы:

 

open #1,"appr.txt","w"

dim v(100),x(100),y(100),f(100)

z=.001

y(1)=0.02:y(2)=0.07:y(3)=0.17:y(4)=0.36:y(5)=0.55:y(6)=0.75

y(7)=0.96:y(8)=1.11:y(9)=1.21:y(10)=1.29:y(11)=1.28:y(12)=1.24

y(13)=1.18:y(14)=1.07:y(15)=0.95:y(16)=0.85:y(17)=0.72:y(18)=0.59

y(19)=0.50:y(20)=0.4:y(21)=0.31:y(22)=0.26:y(23)=0.19:y(24)=0.14

y(25)=0.12:y(26)=0.08:y(27)=0.05:y(28)=0.05:y(29)=0.03:y(30)=0.02

y(31)=0.02:y(32)=0.01

for i=1 to 32:x(i)=i/5:

print i,x(i),y(i)

next i

a0=1:b0=1:c0=1

s1=10^100:nn=10000000

for j=1 to nn

a=a0*(1+z*(ran()-.5))

b=b0*(1+z*(ran()-.5))

c=c0*(1+z*(ran()-.5))

s=0

for i=1 to 32

x=x(i)

f(i)=x^a*b^(- x^c)

s=s+(y(i)-f(i))^2

next i

if s<=s1 then

print a,b,c,s

print #1,a using "#####.########",b using "#####.########",

c using "#####.########",s using ”#####.#########"

s1=s

a0=a:b0=b:c0=c

fi

next j

 

Ну, кажется, нигде не ошибся. Основная идея этой программы состоит в следующем. Как и в методе итераций Ньютона, нужно задаться начальными значениями параметров. Для них находим сумму всех квадратичных отклонений от экспериментальных точек  S. Далее случайным образом изменяем параметры  (эти команды выделены красным цветом). Если новая случайная комбинация параметров дает  S  меньшее, чем предыдущее, то эти параметры оставляем и снова их случайно изменяем. Если же улучшения нет, то вариант изменений пропускаем. Удобство данного метода в том, что можно рассматривать аппроксимирующие функции любой сложности, а количество независимых параметров  ничем не ограничено. Разве только временем работы компьютера. Когда я писал диссертацию, то мне приходилось рассматривать уравнения с двадцатью независимыми параметрами!

- Я тоже буду писать диссертацию!

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

 

Первый столбец – это параметр  a , второй столбец – параметр  b  и третий столбец – параметр  c.  Последний, четвертый столбец – сумма квадратичных отклонений  S . Как видишь,  последняя  S=0.0017878… достаточно мала и можно смело говорить о хорошей сходимости результатов.

Окончательная аппроксимация выглядит так:

 

 

 

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

 

with(plots):X:=[.2,.4,.6,.8,1,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6,2.8,3,3.2,3.4,3.6,3.8,4,4.2,4.4,4.6,4.8,5,5.2,5.4,5.6,5.8,6,6.2,6.4]:Y:=[.02,.07,.17,.36,.55,.75,.96,1.11,1.21,1.29,1.28,1.24,1.18,1.07,.95,.85,.72,.59,.5,.4,.31,.26,.19,.14,.12,.08,.05,.05,.03,.02,.02,.01]:g:=plot([[X[i],Y[i]]$i=1..32],x=0..6.5,style=POINT,symbol=CIRCLE,color=blue):g1:=plot(x^(2.78961755  )*1.81546788^(-x^(1.5021827)),x=0..6.5,y=0..1.4,thickness=3):display(g1,g);

 

Что скажешь, Андрюшенька?  Видишь, как нам удалось почти идеально подобрать кривую, то есть осуществить аппроксимацию?

- Да, получилось красиво! Но  только  формула похожа на паровоз.

- Эх, дорогой мой внучек! В математике часто важно не находить ответ в простом виде, а вообще обнаружить хоть какой-нибудь  безупречный ответ. Запомни это правило на всю жизнь.

- Запомню, дедуля. Аппроксимацию тоже запомню на всю жизнь.

- Знаешь что, Андрюша,  хорошо бы нам попробовать взять какие-нибудь данные из популярной книги и попробовать аппроксимировать лучше, чем в этой книге.  Я сейчас возьму талмуд  Владимира Дьяконова «Maple 8 в математике, физике и образовании» и попытаюсь найти пример. Ну вот, первый попавший рисунок  5.28. на странице 357:

 

 

with(plots):X:=[1,3,4.2,5,6,7,8,9.3,10]:Y:=[2,5,6.2,9,13,16,23,34,42]:g:=plot([[X[i],Y[i]]$i=1..9],x=0..11,style=POINT,symbol=CIRCLE,color=blue):g1:=plot(1.638921889*exp(x*0.329253712755816884),x=0..11,y=0..42,thickness=3):display(g1,g);

 

 

Сумма квадратичных отклонений  S= 7.892687190

 

Как улучшить этот результат? Логично поступить так: корректировать структуру аппроксимирующей функции. Дело это тонкое, требует большого внимания и оригинального мышления. Надо чувствовать формулу, как скульптор чувствует глину.  Вот, Андрюша, кажется нащупал оптимальную структуру. Сейчас слегка подкорректирую программу на языке Yabasic:

 

open #1,"djak.txt","w"

dim v(100),x(100),y(100),f(100)

z=.00001

y(1)=2:y(2)=5:y(3)=6.2:y(4)=9:y(5)=13:y(6)=16:y(7)=23:y(8)=34:y(9)=42

x(1)=1:x(2)=3:x(3)=4.2:x(4)=5:x(5)=6:x(6)=7:x(7)=8:x(8)=9.3:x(9)=10

a0=1:b0=1:c0=1:k0=1

s1=10^100:nn=10000000

for j=1 to nn

a=a0*(1+z*(ran()-.5))

b=b0*(1+z*(ran()-.5))

c=c0*(1+z*(ran()-.5))

k=k0*(1+z*(ran()-.5))

s=0

for i=1 to 9

x=x(i)

f(i)=a+x^b*c^(x^k)

s=s+(y(i)-f(i))^2

next i

if s<=s1 then

print a,b,c,k,s

print #1,a using "#####.########",b using "#####.########",c using "#####.########",k using "#####.########",s using "#####.#########"

s1=s

a0=a:b0=b:c0=c:k0=k

fi

next j

 

Поиск четырех параметров уравнения происходит, слава Богу, при стремительном уменьшении значения S.  А вот и уперлись в экстремум. Андрюша! Набивай очередную графику.

- Говори, дедуля.

 

with(plots):X:=[1,3,4.2,5,6,7,8,9.3,10]:Y:=[2,5,6.2,9,13,16,23,34,42]:g:=plot([[X[i],Y[i]]$i=1..9],x=0..11,style=POINT,symbol=CIRCLE,color=blue):g1:=plot(0.93115383+x^(0.73443836)*1.14406199^(x^(1.17744084)),x=0..11,y=0..42,thickness=3):display(g1,g);

 

- Здесь, в командах системы Maple, уже включена  аппроксимирующая формула. Ее общий вид давай выведем на рисунок:

 

Полученная кривая идеально сгладила все неровности экспериментальных точек. Сумма квадратичных отклонения оказалась рекордно малой:  S = 1.797081329, что в 4 с лишним раза меньше, чем в образцовом книжном решении!

Ну как, Андрюша, понравилось тебе исследование?

- Понравилось! Наша красная кривая, словно змея, ползет между черными точками. В книге не так! В книге точки больше разгуливают туда-сюда.

- Именно в этом и сила метода Монте-Карло. В программе можно экспериментировать, принимая какую угодно структуру  и сколько угодно коэффициентов.

- Все, деда. Пойдем смотреть мультики и играть в лото.

 

10  февраля 2011 г.