Статья Автор: Лебедев Дмитрий

Нахождение лучшего решения уравнения в рациональных числах.

Постановка задачи.
Непрерывная функция f(x) на цлелевом интервале (L,R) возрастает и f(L) <0 < f(R)
Ноебходимо найти корень уравнения f(x) = 0.
В качестве ответа вывести наилучшее "левое приблежение" рациональной дробью со знаменателем не превосходящим N
Пара чисел a, b (a - целое, b - натуральное) будет ответом, тогда и только тогда, когда:

\((I)\ \ L < \frac{a}{b}<R ;\ b<=N;\ f(\frac{a}{b})<0; a\in\mathbb{Z};\ b\in\mathbb{N};\ НОД(|a|,b)=1\)
\(для\ любых\ c,d\ (\frac{c}{d} \neq \frac{a}{b}\ ),\ удовлетворяющих\ (I)\ следует,\ что\ \frac{c}{d}\ <\ \frac{a}{b} \)

Рассмотрим простой пример: \(f(x)=x-\sqrt{m}; N=10^k\)

1, Попробуем положить x=m**0.5 и "округлим знаменатель дроби до N"
Запустив программу для разных m, обнаружим "стабилизацию" точности и округления слева и справа.

from fractions import Fraction as Fr# импорт из модуля рациональных чисел
m=int(input()) # значение для оценки корня
x=Fr(m**0.5) # перевод вещественного в дробь
print(x) # вывод дроби
N,k=1,0 # начальные для округления
while k<20 : # шаги округления
  N,k=N*10,k+1 # значения для округления
  y=x.limit_denominator(N) # результат округления
  print('k=',k,y,float(y*y-m)) # вывод округления

2. Для нахождения корня возмем бинарный поиск на 100 шагов и определим "левую" и "правую" границы и сравним с результатом простого вычисления.
Если посмотреть на результаты для m=2 (корня из двух), то можно увидеть, что совпадение приближений левой и правой границ.
Возникает вопрос: Как найти нужное приблежение?

# Бинарный поиск для корня из m
from fractions import Fraction as Fr # импорт из модуля рациональных чисел
def f_frac(M,C): # вычисление знака целевой функции
    y=M*M-C # значение целевой функции
    if y>0 : return 1 # случай знака +
    if y<0 : return -1   # минус
    return 0 # точное значение
def bin_frac(C,k): # Бинарный поиск для дробей
    C=Fr(C,1) # параметр в дробь
    L,R=Fr(0,1),C # задание целевого интервала
    p=Fr(1,2) # 0,5
    while k>0 : # выполнение шагов алгоритма
        k-=1 # изменение счётчика
        M=(L+R)*p # вычисление середины
        y=f_frac(M,C) # значение целевой функции
        if y==0 : return M,M
        if y>0 :# изменение границ целевого интервала по знаку
            R=M
        else: 
            L=M
    return L,R #  возврат границ целевого отрезка
# Основная программа
m=int(input()) # значение для оценки корня
x=Fr(m**0.5) # перевод вещественного в дробь
xL,xR=bin_frac(m,100) # значения границ целевого интервала 
N,k=1,0 # начальные для округления
while k<20 : # шаги округления
  N,k=N*10,k+1 # значения для округления
  y=x.limit_denominator(N) # результат округления
  yL=xL.limit_denominator(N) # результат округления
  if yL*yL<m : sLR='L' # c "какой стороны" от корня
  else :sLR='R'
  yR=xR.limit_denominator(N) # результат округления
  if yR*yR>m : sLR+='R' # c "какой стороны" от корня
  else :sLR+='L'
  print('k=',k,sLR,'yL=',yL,'y=',y,'yR=',yR) # вывод округления

3. Попробуем изменить способ получения среднего значения на метод "средней скорости".
Поиск продолжается до достижения требуемого значения знаменателя. Эта дробь и будет искомой с нужной границей.
Такой способ подходит для нахождения корней.

# "Средняя скорость" для корня из m
from fractions import Fraction as Fr # импорт из модуля рациональных чисел
def f_frac(M,C): # вычисление знака целевой функции
    y=M*M-C # значение целевой функции
    if y>0 : return 1 # случай знака +
    if y<0 : return -1   # минус
    return 0 # точное значение
def bin_frac(C,k): # Бинарный поиск для дробей
    C=Fr(C,1) # параметр в дробь
    L,R=Fr(0,1),C # задание целевого интервала
    p=Fr(1,2) # 0,5
    while k>0 : # выполнение шагов алгоритма
        k-=1 # изменение счётчика
        M=(L+R)*p # вычисление середины
        y=f_frac(M,C) # значение целевой функции
        if y==0 : return M,M
        if y>0 :# изменение границ целевого интервала по знаку
            R=M
        else: 
            L=M
    return L,R #  возврат границ целевого отрезка
def bin_fracM(C,N): # Бинарный поиск для дробей по Медианту 
    C=Fr(C,1) # параметр в дробь
    L,R=Fr(0,1),C # задание целевого интервала
    LL,RR=[],[]
    k=1
    while 1 : # выполнение шагов алгоритма
        k+=1 # изменение счётчика
        # среднее по правилу (a+c)/(b+d)
        M=Fr(L.numerator+R.numerator,L.denominator+R.denominator)
        if M.denominator > N : # прекращаем при достижении точности
          return LL,RR,k
        y=f_frac(M,C) # значение целевой функции
        if y==0 :
          return [M],[M],k
        if y>0 :# изменение границ целевого интервала по знаку
          R=M
          RR.append(R)
        else: 
          L=M
          LL.append(L)

# Основная программа
m=int(input()) # значение для оценки корня
x=Fr(m**0.5) # перевод вещественного в дробь
kk=20 #  максимальная степень точности
xL,xR=bin_frac(m,100) # значения границ целевого интервала
LL,RR,kz=bin_fracM(m,10**kk) #  получение оценок со "средней скоростью"
N,k=1,0 # начальные для округления
iz,jz=0,0 # начальные индексы z-решения
while k<kk : # шаги округления
  N,k=N*10,k+1 # значения для округления
  while (iz<len(LL)-1) and (LL[iz+1].denominator<=N) : iz+=1
  while (jz<len(RR)-1) and (RR[jz+1].denominator<=N) : jz+=1
  zL,zR=LL[iz],RR[jz]
  y=x.limit_denominator(N) # результат округления
  yL=xL.limit_denominator(N) # результат округления
  if f_frac(yL,m) <= 0 : sLR='L'
  else :sLR='R'
  yR=xR.limit_denominator(N) # результат округления
  if f_frac(yR,m) >= 0 : sLR+='R'
  else :sLR+='L'
  if sLR=='LR' : print('k=',k,sLR,'zL=',zL,'zR=',zR,float(zL-yL),float(yR-zR)) # вывод округления

4. Попробуем решить уравнение \(x^2+\sqrt{x}\ =\ C\)
Возьмем решение и изменим целевую функцию.
Для получения оптимальных решений желательно, чтобы длина целевого интервала была равна 1.
Для этого достаточно получить прикидку бинарным поиском
Объяснение этому можно найти в теории цепных дробей

# "Средняя скорость" для уравнения x*x+sqrt(x)=C
from fractions import Fraction as Fr # импорт из модуля рациональных чисел
def f_frac(M,C): # вычисление знака целевой функции
    y=C-M*M # получение знака
    if y<0 :return 1 # случай знака +
    y=y*y-M
    if y<0 : return 1 # случай знака +
    if y>0 : return -1 # случай знака 1
    return 0 # нуль 
def bin_frac(C,k): # Бинарный поиск для дробей
    C=Fr(C,1) # параметр в дробь
    L,R=Fr(0,1),C # задание целевого интервала
    p=Fr(1,2) # 0,5
    while k>0 : # выполнение шагов алгоритма
        k-=1 # изменение счётчика
        M=(L+R)*p # вычисление середины
        y=f_frac(M,C) # значение целевой функции
        if y==0 : return M,M
        if y>0 :# изменение границ целевого интервала по знаку
            R=M
        else: 
            L=M
    return L,R #  возврат границ целевого отрезка
def bin_fracM(C,N,a=0): # Бинарный поиск для дробей по Медианту 
    C=Fr(C,1) # параметр в дробь
    L,R=Fr(0,1),C # задание целевого интервала
    if a>0 : L,R=Fr(a,1),Fr(a+1,1) # задание целевого интервала
    LL,RR=[L],[R]
    k=1
    while 1 : # выполнение шагов алгоритма
        k+=1 # изменение счётчика
        # среднее по правилу (a+c)/(b+d)
        M=Fr(L.numerator+R.numerator,L.denominator+R.denominator)
        if M.denominator > N : # прекращаем при достижении точности
          return LL,RR,k
        y=f_frac(M,C) # значение целевой функции
        if y==0 :
          LL.append(M)
          RR.append(M)
          return [M],[M],k
        if y>0 :# изменение границ целевого интервала по знаку
          R=M
          RR.append(R)
        else: 
          L=M
          LL.append(L)

# Основная программа
m=int(input()) # значение для оценки корня
x=Fr(m**0.5) # перевод вещественного в дробь
kk=20 #  максимальная степень точности
xL,xR=bin_frac(m,100) # значения границ целевого интервала
print(float(xL))
LL,RR,kz=bin_fracM(m,10**kk,int(xL)) #  получение оценок со "средней скоростью"
N,k=1,0 # начальные для округления
iz,jz=0,0 # начальные индексы z-решения
while k<kk : # шаги округления
  N,k=N*10,k+1 # значения для округления
  while (iz<len(LL)-1) and (LL[iz+1].denominator<=N) : iz+=1
  while (jz<len(RR)-1) and (RR[jz+1].denominator<=N) : jz+=1
  zL,zR=LL[iz],RR[jz]
  y=x.limit_denominator(N) # результат округления
  yL=xL.limit_denominator(N) # результат округления
  if f_frac(yL,m) <= 0 : sLR='L'
  else :sLR='R'
  yR=xR.limit_denominator(N) # результат округления
  if f_frac(yR,m)>=0 : sLR+='R'
  else :sLR+='L'
  if sLR=='LR' : 
    print('k=',k,sLR,'zL=',zL,'zR=',zR,float(zL-yL),float(yR-zR)) # вывод округления

Чтобы оставить комментарий нужна авторизация

Печать