Аннотация
BC - это программа-калькулятор командной строки и простейший язык программирования
арифметических вычислений с цифрами произвольной точности. Возможны манипуляции с числами
длиной более чем несколько тысяч цифр. На языке BC могут быть определены и сохранены
числовые функции для последующего их выполнения.
Содержание
Список примеров
- Преобразование треугольник-звезда
- Расчёт платежей по кредиту
- Индуктивность катушки
- Процентная ставка
- Скорость истечения газа под давлением
- Точка Фейнмана
- ln(2) в двоичной системе
- Округление
- Алгоритм Евклида
- Аннуитетные платежи
- Расчёт вращательной статсуммы молекулы HD
- Расчёт термо-э.д.с. спая Ni-Cr/Ni
- Рекурсивный алгоритм Евклида
- Установка корректного контекста
- Комплексные числа
Введение
Большинство расчётов, которые встречаются в жизни, не требуют продвинутое
программное обеспечение. Использовать для них электронные таблицы MS Excel,
Origin и другие пакеты - всё что палить из пушки по воробьям. Тут достаточно
ручки, бумаги и обычного калькулятора. С другой стороны, а зачем ручка и
бумага? Может быть лучше решать задачу сразу в калькуляторе?
Современные научные калькуляторы уже освоили
Super Visual Perfect Algebraic Method
т.е. записывают свои вычисления формульной строкой. Следовательно, программе на
компьютере нет никакого смысла притворяться старым добрым кнопочным
устройством.
Калькулятор в режиме команд-формул это ничто иное, как
язык программирования.
Пугаться не стоит, всё на самом деле просто. Даже
BASIC,
считающийся языком программирования для начинающих, и тот гораздо сложнее, чем BC.
Язык BC имеет
структуры управления
и режим немедленного выполнения операций. Могут быть определены и сохранены функции для последующего
их выполнения.
Программируемость позволяет BC решать задачи, требующие выбор по условиям, итерации, рекурсии и
другие приёмы, недоступные простым формулам. BC не универсальный язык программирования широкого
применения, как C или C++, однако, это не минус, а скорее плюс. Для своей области – рутинные
инженерные и научные расчёты он очень хорош, а потеря универсальности с лихвой компенсируется
простотой изучения языка.
Калькулятор BC обычно входит в комплект к операционным системам Linux и Mac OS X,
его запуск осуществляется из терминала.
username@computer:~$ bc
Для ОС Windows нужна установка. Скачайте
setup.exe и
устанавливайте. Для быстрого доступа создайте ярлык на программу (обычно она в папке
...\GnuWin32\bin\bc.exe).
Запуск: расчёты по формулам
Пример 1.1. Преобразование треугольник-звезда
Преобразование треугольник-звезда
позволяет упростить расчёт электрических цепей, содержащих замкнутые контуры из резисторов.
Преобразование осуществляется по формулам:
$$
\begin{equation}
\begin{split}
R_1& =\frac{R_{12}R_{13}}{R_{12}+R_{13}+R_{23}}\\
R_2& =\frac{R_{12}R_{23}}{R_{12}+R_{13}+R_{23}}\\
R_3& =\frac{R_{23}R_{13}}{R_{12}+R_{13}+R_{23}}
\end{split}
\end{equation}
$$
Полезный рецепт? Запишем процедуру расчёта в текстовом файле и будём использовать!
Запустим bc с ключом
-l, введём команды (текст внутри
/* и
*/ –
это комментарии) и произведём вычисления.
username@computer:~$ bc -l
/* Входные данные: сопротивления, Ом. */
r12=2; r13=3; r23=5;
rsum=r12+r13+r23
r1=r12*r13/rsum
r2=r12*r23/rsum
r3=r23*r13/rsum
/* Ответ: */
r1
.60000000000000000000
r2
1.00000000000000000000
r3
1.50000000000000000000
Пример 1.2. Расчёт платежей по кредиту
Стиральная машина стоимостью 30000 рублей куплена в кредит под 20% на 12 месяцев.
Погашение долга происходит дифференцированными платежами - ежемесячный платёж
уменьшается со временем, состоит из постоянного долга (неизменный размер) и
процентов на невыплаченный остаток займа. Ежемесячно будет возвращаться 1/12 часть
кредита + начисленные проценты. Задача - рассчитать график платежей.
/* Обозначения:
sum - сумма кредита,
pay - ежем. платёж,
int - процентная ставка,
rem - остаток по кредиту
*/
sum=30000; int=0.20;
/* первый платёж */
rem=sum;
/* 1/12 кредита + % от долга */
pay=sum/12 + rem*int/12; pay
3000.00000000000000000000
rem=rem-sum/12; rem
27500.00000000000000000000
/* следующий платёж */
pay=sum/12 + rem*int/12; pay
2958.33333333333333333333
rem=rem-sum/12; rem
25000.00000000000000000000
/* и так далее */
pay=sum/12 + rem*int/12; pay
2916.66666666666666666666
rem=rem-sum/12; rem
22500.00000000000000000000
Пример 1.3. Индуктивность катушки
Требуется рассчитать индуктивность однослойной катушки квадратного сечения, см. рисунок ниже.
Находим индуктивность по формулам:
$$
\begin{equation}
\begin{split}
F& =2h\ln \frac{D_{2}}{D_{1}}\\
L& =F\omega^2
\end{split}
\end{equation}
$$
/*
Входные данные: высота катушки,
внутренний и наружный диаметры,
число витков.
*/
h=1; d1=1; d2=3; omega=100;
/* Формула расчёта индуктивности */
f=2*h*l(d2/1)
f*omega^2
21972.24577336219382780000
Убедитесь в том, что вызываете BC с ключом
-l, т.е.
bc.exe -l.
Ключ
-l выставляет точность по умолчанию
scale=20 цифр
после запятой и подключает функции математической библиотеки:
s(x)- $\sin x$, аргумент $x$ в радианах
c(x)- $\cos x$, аргумент $x$ в радианах
a(x)- $\arctan x$, результат вычисления в радианах
l(x)- $\ln x$ натуральный логарифм
e(x)- $e^x$ экспонента
j(n, x)- $J_n(x)$ - функции Бесселя
Без ключа
-l в BC доступны только совсем элементарные действия:
сложение +, вычитание -, умножение *, деление /, нахождение остатка %, возведение в
целочисленную степень ^ и извлечение квадратного корня
sqrt(x).
По умолчанию
scale=0, потому
sqrt(2) вернёт 1, а не 1.414213562…
Следовательно, более-менее привычный режим калькулятора начинается с команды
scale=n,
где n – требуемая точность. Однако, остаток от деления a на b ищется по формуле
$$
\begin{equation}
a\% b=a-(a/b)\cdot b
\end{equation}
$$
следовательно, правильный результат в этом случае только при
scale=0.
Функции и переменные
Пример 2.1. Процентная ставка.
Необходимо привести процентную ставку 20% из годового значения к месячному, т.е. вычислить
$\sqrt[12]{100\%+20\%}-1$.
Для решения нам понадобится операция $y^x$, которой BC не располагает.
Хорошо, определим нужную нам функцию из имеющихся:
$$
\begin{equation}
y^x = e^{y\cdot\ln x}
\end{equation}
$$
На языке BC это делается по следующему шаблону.
/*
define функция(аргумент1, аргумент2 … ) {
return результат
}
*/
define pow(x,y) {
return e(l(x)*y)
}
pow(1.20,1/12)
.01530947049973121704
Пример 2.2. Скорость истечения газа под давлением.
Рассмотрим формулу для вычисления скорости истечения идеального газа из сосуда под давлением.
Скорость истечения зависит от давления внутри и вне сосуда, от молекулярной массы газа, от
температуры и (из-за адиабатического расширения) от отношения молярных теплоёмкостей
$k=C_p/C_v$. Формула имеет вид
$$
\begin{equation}
v=\sqrt{\frac{2RT}{\gamma M}\left[1-\left(\frac{p}{p_{0}}\right)^{\gamma}\right]}
\end{equation}
$$
где $v$ – скорость истечения; $R$ – универсальная газовая постоянная; $T$ – температура
в сосуде; $p$ – внешнее давление; $p_0$ – давление в сосуде; $M$ – молекулярная масса газа;
$\gamma$ – вспомогательная величина, равная $\gamma=(k-1)/k$.
/* Температура в К */
t=293
/* Молекулярная масса, кг/моль */
m=28.96*10^-3
/* Отношение k=Cp/Cv */
k=1.402
/* Внешнее давление */
p=1.2
/* Внутреннее давление */
p0=2.4
/* Универсальная газовая постоянная, Дж/К/моль */
r=8.31441
g=(k-1)/k
/* Формула, м/с */
v=sqrt(2/g*r*t/m*(1-pow(p/p0,g)))
/* Ответ */
v
325.20441311481813026525
Подытожим.
- Имена переменных
- В именах переменных только латинские буквы, и только прописные.
Выражения
A=1 или длина=10.0 не допускаются правилами.
- Комментарии
- Многострочные комментарии должны быть заключены в
/* ... */,
а короткие, длиной в одну строку должны начинаться с символа #
- Десятичный разделитель
- Только точка. Никаких 3,14.
Пример 2.3. Точка Фейнмана.
Точка Фейнмана – последовательность
из шести девяток, начинающаяся с 762-ой цифры десятичной записи числа $\pi$.
Носит имя американского физика Ричарда Фейнмана, который сказал на одной лекции,
что хотел бы запомнить цифры числа $\pi$ до этой позиции, чтобы заканчивать рассказ
кому-либо словами «девять, девять, девять, девять, девять, девять и так далее»,
как бы предполагая, что значение $\pi$ рационально. Проверим это!
/* BC способен работать с числами произвольной точности */
scale=768
4*a(1)
3.1415926535897932384626433832795028841971693993751058209\
749445923078164062862089986280348253421170679821480865132\
823066470938446095505822317253594081284811174502841027019\
385211055596446229489549303819644288109756659334461284756\
482337867831652712019091456485669234603486104543266482133\
936072602491412737245870066063155881748815209209628292540\
917153643678925903600113305305488204665213841469519415116\
094330572703657595919530921861173819326117931051185480744\
623799627495673518857527248912279381830119491298336733624\
406566430860213949463952247371907021798609437027705392171\
762931767523846748184676694051320005681271452635608277857\
713427577896091736371787214684409012249534301465495853710\
507922796892589235420199561121290219608640344181598136297\
74771309960518707211349999996
Пример 2.4. ln(2) в двоичной системе.
Алгоритм «текущего крана» (англ. –
Spigot Algorithm)
позволяет находить числа после запятой в константах, таких как $\pi$, $e$, $\ln(2)$ и других,
не вычисляя предыдущие цифры. То есть, если нужно найти 762-ую цифру в десятичной записи
$\pi$ – алгоритм найдёт её не вычисляя предыдущие цифры с самой первой и по 762.
Например, с помощью соотношения
$$
\begin{equation}
\ln(2)=\sum_{k=1}^{\infty}\frac{1}{k2^{k}}
\end{equation}
$$
удаётся находить цифры ln(2) в двоичной записи начиная с n-той позиции.
В статье Wikipedia
Spigot Algorithm
утверждается, что
$$
\begin{equation}
2^{7}\ln(2)\mod1\approx\frac{64}{105}+\frac{37}{360}\approx0.1011\ldots_{2}
\end{equation}
$$
и мы проверим это с помощью BC.
l(2)*2^7
88.7228391040
obase=2
last
1011000.1011100100001011111110111100011111
/* цифры после запятой: 1011... */
64/105+37/360
.1011011001011001011001011001010111
obase=A
l(2)*2^7
88.7228391040
Подытожим.
- Переменная
obase
- (output base) устанавливает какую систему использовать при отображении
результатов,
obase=2 приводит к тому, что все результаты вычислений
выводятся в двоичной системе. Формально корректное значение obase=1
BC считает ошибочным.
- Переменная
ibase
- (input base) устанавливает основание
для ввода чисел.
- Переменная
last
- хранит последний вычисленный результат.
- Символы верхнего регистра: A, B, C, D…
- обозначают числа 10, 11, 12,
13… в системах выше десятичной, независимо от установок obase/ibase, A это всегда 10,
что позволяет командой
obase=A вернуться назад в десятичную систему
Алгоритмы
Логические операции
Допустим, мы хотим задать функцию вычисления модуля числа $f(x)=|x|$.
Если x>0 то возвращаем x, а если нет – то возвращаем -x. Как это сделать в BC?
Рассмотрим всё по-порядку.
Сравнение чисел –
логическая операция, как и сложение-вычитание
арифметическая операция. Истинное утверждение возвращает 1, а ложное – 0.
x < y- $x<y$
x <= y- $x\leq y$
x > y- $x>y$
x >= y- $x\geq y$
x == y- $x=y$
x != y- $x\neq y$
Заметим, что
присваивание = и
сравнение на равенство ==
обозначены по-разному, так как BC не видит контекст употребления символа
=.
В общем случае, если мы хотим указать диапазон допустимых значений переменной x,
нам нужны не только логические утверждения (сравнения), но и операции, такие как И, ИЛИ, НЕ.
Они обозначаются в BC знаками
&&,
|| и
!,
соответственно.
Скажем, утверждение $x\in\left(a,b\right]$ на языке BC:
x > a && x<= b
Отрицание этого утверждения записывается как
!(x > a && x<= b)
Приоритеты операций
В последнем примере логического утверждения, что
$x$ не содержится в интервале
$\left(a,b\right]$ нам понадобилось заключить выражение
x > a && x<= b
в скобки. В самом деле, иначе было бы непонятно, отрицание какого именно утверждения производит оператор
!.
Можно было бы предположить, что выражение
!x > a && x<= b
означает отрицание первого утверждения,
x > a и операцию логического И со
вторым утверждением
x<= b. Однако, это не так.
Обратите внимание, что истинное утверждение равно 1, а ложное – 0, однако,
истине соответствует также любое ненулевое значение.
!-0.1
0
!100
0
!0
1
Как видно, логическое НЕ от любой величины не равной нулю всегда 0, то есть ложь.
Так что !x > a означает не отрицание утверждения x > a,
а сравнение отрицания x с числом a.
Такие у BC правила! Иначе пришлось бы вводить разделение переменных на типы –
булевый
(логический, 0 или 1) и арифметический, что многое бы усложнило.
Язык BC безтиповый (о том, что это значит в языках программирования читайте
здесь).
Чтобы не запутаться, нужно ввести
приоритеты операций. Рассмотрим их в порядке от высшего к низшему приоритету.
- Инкремент/декремент
++ и -- (они увеличивают/уменьшают значение переменной на 1, о них подробнее позже). Для x=0.1 x++ вернёт 1.1 и увеличит значение x на единицу;
- Оператор
- перед числом или переменной. Если x=2, то выражение -x вернёт -2;
- Возведение в степень
^;
- Умножение, деление и остаток по модулю
*,/, %;
- Сложение и вычитание
+ и -;
- Присваивание
=;
- Сравнения
>,<, >=, <=, == и !=;
- Отрицание (логическое НЕ)
!;
- Логическое И
&&;
- Логическое ИЛИ
|| .
Согласно правилам математики, $-x^2<0$, тогда как в языке BC (и внезапно, в MS Excel)
это не так. Оператор
- с одним аргументом ( как -x, -1.0)
имеет более высокий приоритет, чем возведение в степень
^. Однако,
простое вычитание с двумя аргументами (как a-b, a-1.0) имеет приоритет ниже,
чем степень.
/* (-3)^2 */
-3^2
9
/* 0-(2)^2 */
0-2^2
-4
Выражение
a = 3 < 5 сделает
a=3, затем сравнит результат
(результат присваивания
x=y всегда
y) с 5 и вернёт 1 (истина).
А как мы узнали, что результат операции присваивания такой?
С помощью специальной переменной
last, которая хранит результат последней операции.
a = 3
last
3
a = 3 < 5
1
Теперь у нас есть возможность вычислить $|x|$:
x*=-1^((x>0)+1)
- Операторы
+=, -=, *=, /=, %=, ^=
- Всё это варианты присваивания типа
x = x + y или x = x^y,
записанные в сокращенной записи x+=y и x^=y,
соответственно.
Как вы думаете, какое значение приобретёт переменная
x в следующем выражении?
x > 0 || x*=-1
Штука вся в том, что это выражение сработает как $|x|$!
BC начинает вычислять сравнение
x > 0, если x – положительное число,
то результат всего выражения
x > 0 || x*=-1 очевиден,
т.к.
истина ИЛИ
истина/ложь равно
истина!
Это свойство функции ИЛИ. Следовательно, можно не вычислять следующий аргумент
логического оператора ИЛИ. Пойти по
короткому пути вычислений,
так сказать.
Управляющие конструкции
Теорема Бёма-Якопини
утверждает, что любой исполняемый алгоритм может быть преобразован к структурированному виду,
то есть такому виду, когда ход его выполнения определяется только при помощи трёх структур
управления: последовательности инструкций, ветвлений и циклов.
Последовательности инструкций
Последовательность инструкций в BC это просто список формул.
/* Входные данные: сопротивления, Ом. */
r12=2; r13=3; r23=5;
rsum=r12+r13+r23
r1=r12*r13/rsum
r2=r12*r23/rsum
r3=r23*r13/rsum
Точка с запятой
; используется для записи последовательности в строчку.
Отдельные ветви последовательностей вычислений, например, в функциях, заключены в фигурные
скобки
{}.
Ветвления if … else
Эта конструкция разветвляет поток команд управления, добавляет действия на случай
предусмотренный условием.
Пример 3.1. Округление.
define round(x,n) {
/*
Округление x до n-того знака.
*/
save = scale
if(x>=0) { sign=1 } else { sign=-1 }
x*=sign
scale=0
/*
Вспомним, что scale должно быть 0,
чтобы вычислить остаток от деления.
Деление на 1 при scale=0 отбросит
дробную часть числа.
*/
digits=x*10^(n+1)/1
ld = digits % 10
if(ld>=5) {
digits+=(10-ld)
/*
Ключевое слово else всегда идёт после }.
Такой вариант ОШИБОЧЕН!
if(ld>=5) {
digits+=(10-ld) }
else {
digits-=ld
}
*/
} else {
digits-=ld
}
scale=n
ans=sign*digits/10^(n+1)/1
scale=save
return ans
}
round(sqrt(2),5)
1.41421
Цикл while
Конструкция
while(условие) инструкции
осуществляет повторение инструкций до тех пор, пока остаётся в силе условие.
Вычисления идут по циклу. Выйти из него можно специальной командой
break.
Перейти к следующей итерации, не завершая предыдущую -
continue.
Пример 3.2. Алгоритм Евклида.
Алгоритм Евклида
позволяет найти наибольший общий делитель двух целых чисел.
define gcd(u,v) {
save=scale
scale=0
/* Вместо v!=0 достаточно
написать просто v.
*/
while(v) {
t=u; u=v; v=t%v
}
scale=save
return u
}
gcd(1071,462)
21
Цикл for
При составлении программ часто повторяется следующий шаблон:
i=1
while(i<n) {
/* операции */
i+=1
}
Для его сокращенной записи в языке BC существует конструкция
for(i=1; i<n; i++){/* операции */}.
Внутри скобок имеются три выражения, разделяемые точкой с запятой. Первое выражение -
инициализация
i=1 выполняется один раз перед тем, как войти в цикл.
Второе - проверка условия продолжения цикла
i<n.
Условие вычисляется, и если оно истинно, выполняется тело цикла.
Затем осуществляется приращение счётчика
i++ и условие вычисляется снова.
Цикл заканчивается, когда условие становится ложным.
Заметим, что цикл
for> - более общая конструкция, чем цикл
while). Так, цикл
while(i<n){...} эквивалентен
for(;i<n;) {...}. А какую известную в комбинаторике функцию
вычисляет этот код?
define fun(n) {
for(f=n; --n; f*=n){}
return f
}
fun(1)
1
fun(5)
120
Пример 3.3. Аннуитетные платежи.
Аннуитетными платежами (то есть равновеликими) называются платежи, равные друг другу
на протяжении всего срока кредита. Структура аннуитетного платежа при этом состоит из
суммы погашения основного долга и процентов за пользование кредитными деньгами.
Соотношение этих величин со временем меняется: увеличивается сумма погашения основного
долга, размер процентов, напротив, увеличивается (в пределах фиксированного аннуитетного
платежа). Важный момент: проценты начисляются всегда на остаток долга, при аннуитетных
платежах вначале сумма основного долга минимальна, а значит, размер процентов заметно
больше и убывает незначительно. В первые месяцы приходится платить в основном проценты по кредиту.
Определить величину аннуитетного платежа $P$ можно
по формуле:
$$
\begin{equation}
P=PV\cdot\frac{r}{1-(1+r)^{-n}}
\end{equation}
$$
где $PV$ - сумма кредита, $r$ - процентная ставка за один месяц, $n$ - количество месяцев
на протяжении всего действия аннуитета.
Рассчитаем выплаты по ипотеке на сумму 1000000 рублей под 20% годовых в течении 5 лет.
Выведем таблицу с колонками: номер платежа, задолженность по кредиту,
начисленные проценты, основной долг.
/* Условия ипотеки */
sum=1000000; r=0.2/12; period=12*5;
/* Формула ежемесячного
аннуитетного платежа.
*/
p=sum*r/(1-(1+r)^-period)
for(n=1; n<=period; n++) {
l=sum*r; d=p-l;
/*
Обратная косая черта - \ символ переноса
строки.
*/
print n,"\t",round(sum,2),"\t", \
round(l,2),"\t",round(d,2),"\n"
sum-=d;
}
В языке BC нет переменных-строк, весь необходимый вывод обеспечивается оператором
print.
Он выводит пояснения, записанные в кавычках. Кроме того, он понимает управляющие символы,
такие как
"\t" – табуляция, перевод строки
"\n" и другие.
Пример 3.4. Расчёт вращательной статсуммы молекулы.
Статистическая сумма (или статсумма) (обозначается $Z$, от нем.
Zustandssumme — сумма по состояниям)
— важная величина в статистической физике, содержащая информацию о термодинамических свойствах системы.
Свободная энергия, энтропия и давление, могут быть выражены через статистическую сумму и её производные
Рассчитаем для справочной таблицы статсуммы молекулы HD при температурах от 10.0 до 300 К, воспользовавшись
справочными данными о моменте инерции молекулы
$I=0.61\cdot10^{-47}\text{кг/м}^2$. Вращательная статсумма равна
$$
\begin{equation}
Z=\sum_{J=0}^{\infty}(2J+1)\cdot\exp(-\frac{J(J+1)\hbar^2}{2IkT})
\end{equation}
$$ где $J$ -- номер уровня, $\hbar\approx1.0545718\cdot10^{-34}\text{Дж}\cdot\text{c}$ -- редуцированная постоянная Планка и $k\approx1.38064852\cdot10^{-23}\text{Дж/K}$ -- константа Больцмана.
Введём характеристическую вращательную температуру HD $$\theta=\frac{\hbar^2}{2Ik}$$ и вычислим статсумму по
формуле
$$
\begin{equation}
Z(T)=\sum_{J=0}^{\infty}(2J+1)\cdot\exp(-J(J+1)\frac{\theta}{T})
\end{equation}
$$
/*
Фундаментальные физические константы
http://physics.nist.gov/constants
k=1.38064852*10^(-23); h=1.054571800*10^(-34);
*/
/* Характеристическая температура */
theta = 1.054571800^2/(2*0.61*1.38064852)*10^(-34*2+47+23)
/* Функция расчёта статсуммы */
define zsum(t,n) {
sum=1.0
for(j=1; j<n; j++) {
p=j*(j+1)*theta/t
if(p>50) break
sum+=(2*j+1)*e(-p)
}
return sum
}
/* Расчёт справочной таблицы */
for(t=10; t<300; t+=10) {
print round(zsum(t,50),8),\
" at T =", t, "\n"
}
К этому моменту мы добрались до составления полноценной программы на языке BC.
Разберём внимательно все "подводные камни".
Первый. Ввод чисел, меньших, чем
10^(-scale)
k=1.38064852*10^(-23); h=1.054571800*10^(-34);
theta = h^2/(2*0.61*10^(-47)*k)
Runtime error (func=(main), adr=30): Divide by zero
Дело в том, что при
scale=20, числа меньшие, чем
10^(-scale)
рассматриваются как нули. А как узнать, какое значение нам нужно?
Функция
scale(0.61*10^(-47)) вернёт необходимую точность.
Второй. Кажется, что без следующей проверки можно обойтись:
p=j*(j+1)*theta/t
if(p>50) break
однако, это не так. Без неё инструкция
sum+=(2*j+1)*e(-p) работает очень долго,
вычисления фактически зависают. Это связано с реализацией функции
e(x) в библиотеке BC.
Пример 3.5. Расчёт термо-э.д.с. спая Ni-Cr/Ni.
Обычно измерение термо-э.д.с. проводят двумя соединёнными навстречу друг другу термопарами, например, термопарами Ni-Cr/Ni. При этом один спай термостатируют при 0C, в то время как другой имеет температуру объекта. Экспериментально полученную зависимость термо-э.д.с. от температуры в определённом температурном интервале обычно описывают полиномом.
Пусть разность э.д.с. холодного и горячего спаев термопары и температура горячего спая связаны следующей зависимостью:
$$
\begin{equation}
\begin{split}
C&=25.4498V-0.559195V^{2}+0.10452439V^{3}-8.776153\cdot10^{-3}V^{4}+\\&+3.76041\cdot10^{-4}V^{5}-8.64943\cdot10^{-6}V^{6}+1.021005\cdot10^{-7}V^{7}-\\&
-4.891009\cdot10^{-10}V^{8}
\end{split}
\end{equation}
$$
где $V$ – термо-э.д.с., мВ; $C$ – температура, $^{\circ}$C.
Задача – определить значение термо-э.д.с. при заданной температуре C=200$^{\circ}$C.
Решение. Используя
метод хорд ,
вычисляем корень уравнения. Для этого выбираем интервал $[x_{i-1},\,x_{i}]$,
содержащий искомый корень, т.е. $f(x_{i-1})<0$ и $f(x_{i})>0$. Затем, через две точки
функции проводим хорду, и её пересечение с осью $x$ даёт приближенное решение уравнения.
Повторяем эту процедуру до сходимости.
$$
\begin{equation}
x_{i+1}=x_{i}-\frac{x_{i}-x_{i-1}}{f(x_{i})-f(x_{i-1})}\cdot f(x_{i})
\end{equation}
$$
define abs(x) { if(x < 0) return -x else return x }
/* Разность э.д.с. и температура горячего спая */
define c(v) {
return 25.4498*v - 0.559195*v^2 + 0.10452439*v^3 \
- 8.7761538*10^-3*v^4 + 3.76041*10^-4*v^5 \
- 8.64943*10^-6*v^6 + 1.021005*10^-7*v^7 \
- 4.891009*10^-10*v^8;
}
/* Подбираем диапазон значений э.д.с. дающих
разумные температуры (-273..1200 С) */
c(-6)
-210.13802929414525440000
c(55)
1111.97762786214843750000
x0=-6; x1=55;
t=200;
/* Ищем корень уравнения методом хорд с ошибкой 0.1% */
for(iter=0; abs(x0-x1)>0.001*abs(x1); iter++) {
f0 = c(x0)-t
f1 = c(x1)-t
xn = x1-f1*(x1-x0)/(f1-f0)
x0 = x1; x1 = xn;
}
/* Затрачено итераций */
iter
4
/* Ответ: */
xn
8.17592012890060407645
/* проверка: */
c(xn)
199.99998979466756124429
Поведение функций и переменных в BC
Область видимости переменных
Рассмотрим код из
примера 3.2
define gcd(u,v) {
save=scale
scale=0
while(v) {
t=u; u=v; v=t%v
}
scale=save
return u
}
Допустим, где-то раньше использовались переменные
save,
t.
После вызова функции
gcd(1071,462) прежние значения переменных потеряются.
А что произойдёт с
u и
v?
save=100; t=-50; u=100500; v=3.14
gcd(1071,462)
save
20
t
147
u
100500
v
3.14
А
u и
v остались прежними!
Переменные, получившие имя как параметры при создании какой-либо функции
(внутри скобок) не переписывают значения
глобальных переменных,
созданных приравниванием.
Функции в BC создают свой список переменных, с чистого листа.
Имя
u внутри функции
gdc() означает что-то другое,
чем
u "снаружи". В свою очередь, для функции используемой
внутри функции
gcd(), переменная
u становится
внешней, но и она может определить свою версию и понимание
u.
Проблема в том, что используемые временно
save и
t не находятся в скобках. В BC внутри функции можно явно указать,
что нужна своя собственная переменная
u, а не внешняя с помощью
ключевого слова
auto которое помещается в первой строчке функции.
define gcd(u,v) {
auto save, t
save=scale
scale=0
while(v) {
t=u; u=v; v=t%v
}
scale=save
return u
}
save=100; t=-50; u=100500; v=3.14
gcd(1071,462)
save
100
t
-50
u
100500
v
3.14
Обратите внимание, что в BC существует два вида переменных –
глобальные
(объявленные приравниванием), своего рода "мировые константы", и
локальные
– впервые появившиеся в параметрах или после ключевого слова
auto.
Что локально на одном уровне, то на следующем уровне вложения (глубже) глобально.
/* функции объявленные со словом
void (англ. - пустой)
не возвращают значений с помощью
return
*/
define void fun() { print t, "\n" }
define gcd(u,v) {
auto save, t
save=scale
scale=0
while(v) {
t=u; u=v; v=t%v
}
fun()
scale=save
return u
}
save=100; t=-50; u=100500; v=3.14
fun()
-50
gcd(1071,462)
147
21
save
100
t
-50
Функции BC под свои параметры берут "чистый лист" и не мешают друг другу. В программировании это называется
передача аргументов по значению.
Следует понимать так, что вызываемой функции посылаются не сами аргументы, а их значения.
Рассмотрим код
x=1; y=2
define void swap(x,y) {
auto t
t=x; x=y; y=t
}
swap(x,y)
x
1
y
2
Ожидалось, что функция
swap(x,y) поменяет местами значения x и y.
Это не произошло. Произошло фактически вот что: была вызвана функция
swap(1,2),
она присвоила своим локальным переменным значения
x=2; y=1 и закончила работу.
Внешние x и y остались без изменений.
Рекурсия
Передача только значений открывает новые возможности функций. Раз каждая имеет дело со своим черновиком,
то она может вызвать саму же себя, передать себе значения и задействовать новые ячейки памяти, не помешав никому.
Пример 4.1. Рекурсивный алгоритм Евклида.
define gcd(u,v) {
auto save
depth = depth + 1
save=scale
if(v) {
return gcd(v, u%v)
} else {
return u
}
scale=save
}
depth=0
gcd(1071,462)
21
depth
4
Рекурсивный алгоритм прошёл
depth=4 уровня вызовов, а значит,
для работы ему потребовалось 8 ячеек памяти: по 4 копии параметров u и v.
Итерационная версия алгоритма Евклида потребовала бы только 3 ячейки памяти.
Это общий недостаток рекурсивных алгоритмов – им требуется больше памяти,
однако, в тех случаях, когда глубина рекурсии небольшая, и они позволяют выразить идею
проще, их стоит применять.
Пример 4.2. Установка корректного контекста.
Мы хотим, чтобы функция всегда работала с десятичными числами.
Рекурсивный вызов решает задачу элегантно.
define decimal(x) {
auto base, value
if (ibase != A) {
base = ibase
ibase = A
value = decimal(x)
ibase = base
return value
}
return x/10
}
ibase=2
decimal(100) /* 100 в двоичной системе это 4 */
.40000000000000000000
Область видимости функций
Функции в BC проектируют дедуктивным методом: от общего к частному. Выделим метод
решения уравнений из
примера 3.5 в отдельную функцию
root(x0,x1,eps). На этом этапе мы не уточняем, какое именно
уравнение будем решать, просто обозначим его
eq(x). Также не
конкретизируем, в каком смысле интервал $(x_n,x_{n+1})$ сжимается, это
позже сделает функция
diff(x0,x1). Именно это поведение BC
называется в программировании
позднее связывание.
define root(x0,x1,eps) {
/*
Функция находит корень уравнения eq(x)=0 методом
https://en.wikipedia.org/wiki/False_position_method
на интервале [x0,x1]. В указаном интервале
функция eq(x) должна пересекать ось x, т.е.
f(x0)<0 и f(x1)>0 ИЛИ f(x0)>0 и f(x1)<0
*/
auto diff,
f0=eq(x0);
f1=eq(x1)
is=((f0 maxit=100 x0="" x1="" xn="">0 && f1<0 f0="" f1="">0))
if(!is) {
print "корень ур-ния не изолирован\n"
return x0
}
for(i=0; diff(x0,x1)>eps; i++) {
f0 = eq(x0)
f1 = eq(x1)
xn = x1-f1*(x1-x0)/(f1-f0)
x0 = x1; x1 = xn;
if(i>maxit) {
print "Cходимость не достигнута за "
print maxit, " итераций\n"
break
}
}
return xn
}
Теперь задача из
примера 3.5 решается так:
define abs(x) {
if(x<0) { return -x
} else { return x }
}
define diff(a,b) { return abs(a-b) }
/* Разность э.д.с. и температура горячего спая */
define c(v) {
return 25.4498*v - 0.559195*v^2 + 0.10452439*v^3 \
- 8.7761538*10^-3*v^4 + 3.76041*10^-4*v^5 \
- 8.64943*10^-6*v^6 + 1.021005*10^-7*v^7 \
- 4.891009*10^-10*v^8;
}
/* Определяем ур-ние, которое предстоит решать
функции root()
*/
define eq(x) {
return c(x)-200
}
/* Подбираем диапазон значений э.д.с. дающих
разумные температуры (-273..1200 С) */
c(-6)
-210.13802929414525440000
c(55)
1111.97762786214843750000
emf=root(-6,55,0.0001); emf
8.17592054055875721094
c(emf)
199.99999999999800301248
Массивы чисел
Калькулятор BC располагает средствами работы с массивами,
т.е. индексированными величинами $x_i$, где $i=0\ldots N$. Воспользуемся
этим для реализации арифметики комплексных чисел.
re=0
im=1
/* запишем комплексное число
a = a0 + a1*i, где i*i=-1
как массив
*/
a[re]=0; a[im]=1;
Далее определим основные функции: сложение, вычитание,
умножение и деление. К сожалению, функции BC не позволяют
возвращать массив.
define add(a[],b[]) {
auto t[]
t[re]=a[re]+b[re]
t[im]=a[im]+b[im]
return t[]
(standard_in) 8: syntax error
}
Следующий вариант тоже не сработает, так как функции BC принимают
только значения аргументов.
re=0; im=1;
define void add(a[],b[],c[]) {
c[re]=a[re]+b[re]
c[im]=a[im]+b[im]
}
a[re]=0; a[im]=1;
b[re]=1; b[im]=0;
add(a[],b[],c[])
c[re]
0
c[im]
0
Плохо, очень плохо. Иногда
очень нужно, чтобы функция могла
изменить аргумент. Вот специально на такие случаи в BC есть лазейка:
передача функции
*a[]. Переменная-массив, помеченная
звёздочкой.
Пример 4.3. Комплексные числа.
re=0; im=1; rad=2; phi=3;
define void set(*a[],b[]) {
a[re]=b[re]; a[im]=b[im]
a[rad]=b[rad]; a[phi]=b[phi]
}
define void add(a[],b[],*c[]) {
c[re] = a[re]+b[re]
c[im] = a[im]+b[im]
}
define void sub(a[],b[],*c[]) {
c[re] = a[re]-b[re]
c[im] = a[im]-b[im]
}
define void mul(a[],b[],*c[]) {
c[re] = a[re]*b[re]-a[im]*b[im]
c[im] = a[im]*b[re]+a[re]*b[im]
}
define void div(a[],b[],*c[]) {
c[re] = a[re]*b[re]+a[im]*b[im]
c[re]/= b[re]^2 + b[im]^2
c[im]=a[im]*b[re]-a[re]*b[im]
c[im]/= b[re]^2 + b[im]^2
}
define void polar(*a[]) {
a[rad] = sqrt(a[re]^2+a[im]^2)
a[phi] = a(a[im]/a[re])
}
define void alg(*a[]) {
a[re] = a[rad]*c(a[phi])
a[im] = a[rad]*s(a[phi])
}
/* Формула Муавра */
define void pow(n,a[],*c[]) {
c[rad] = a[rad]^n
c[phi] = a[phi]*n
alg(c[])
}
define void cprint(a[]) {
if(a[im]>0) {
print a[re]," + i*",a[im],"\n"
} else {
print a[re]," - i*",abs(a[im]),"\n"
}
}
/* создаём комплексное число */
a[re]=1; a[im]=1; polar(a[])
cprint(a[])
1 + i*1
/* проверим формулу Муавра */
define void ipow(n, a[],*c[]) {
set(c[],a[])
while(--n) mul(c[],a[],c[])
}
ipow(5,a[],c[])
cprint(c[])
-4 - i*4
pow(5,a[],c[])
cprint(c[])
-4.00000000000000000008 - i*3.99999999999999999985
Дополнительная информация
Калькулятор нужен от случая к случаю. За это время можно многое забыть.
Тем и хорош BC - не придётся долго вспоминать, он очень простой. На всякий случай,
полезно сделать себе шпаргалку, прямо в BC. Вызов -
help().
define void help() {
print "/* Комментарии */\n"
print "\tАрифметика:\n"
print "+,-,/, sqrt() - корень квадратный\n"
print "% (остаток), ^ - степень x^n = x*...*x\n"
print "\tВызов\n"
print "bc -l, scale=20 (точность)\n"
print "a(x) - арктангенс, s/c - синус/косинус\n"
print "l(x) - логарифм, e(x) - экспонента\n"
print "\tЛогика\n"
print "0 - ложь, 1 - истина\n"
print "==, !=, <=, >=, < > - сравнения\n"
print "!x, x && y, x || y - НЕ, И, ИЛИ\n"
print "\tВвод/вывод\n"
print "ibase/obase - основания чисел для\n"
print "ввода/вывода\n"
print "x=read()/print x - чтение/печать\n"
print "\tФункции\n"
print "define [void] fun(x,y) {\n"
print "auto t, z\n"
print "/* временные переменные - auto */\n"
print "return answer }\n"
print "\tКонструкции алгоритмов\n"
print "Выбор: if(условие) ... else ...\n"
print "Циклы: while(условие) { ... }\n"
print "for(i=0; i<n; i++) { ... }\n"
print "break/continue - выход из цикла/\n"
print "продолжить итерации\n"
print "\tМассивы\n"
print "fun(a[]); a[0]=1.0\n"
print "с изменением аргумента: fun(*a[])\n"
}
Доскональные знания даёт опыт работы и
Руководство пользователя.
Что касается опыта, то профессор математики
Keith Matthews, создал много
примеров использования
BC в современной теории чисел.
Удачи в вычислениях!
Пусть в кубическом сантиметре воды N3 молекул, тогда на поверхности
кубика – N2.
В ходе нарезания кубика исчезло N2 связей между молекулами. Полагая, что с
каждой разорванной связью исчезает энергия Eb, мы получаем, что