Быстрая и грубая оценка количества значащих двоичных цифр в заданном неотрицательном числе

#алгоритм #численные_методы

Собственно, знаю о решении в лоб - двоичном логарифме, и итерационной его реализации
для целого результата: D сдвигов вправо, пока заданное число N > 0:
int D = 0;
while (N > 0) {
    N = N >> 1;
    D++;
}

Но существует ли какой-либо хак, чтобы получить результат быстрее, чем сдвигами в
цикле? (возможно, пожертвовав точностью)    

Ответы

Ответ 1

Количество нулевых бит в 64-разрядном целом
// adapted from Hacker's Delight
int clzll(uint64_t x) {
  int n;

  if (x == 0) return(64);
  n = 0;
  if (x <= 0x00000000FFFFFFFFL) {n = n + 32; x = x << 32;}
  if (x <= 0x0000FFFFFFFFFFFFL) {n = n + 16; x = x << 16;}
  if (x <= 0x00FFFFFFFFFFFFFFL) {n = n + 8; x = x << 8;}
  if (x <= 0x0FFFFFFFFFFFFFFFL) {n = n + 4; x = x << 4;}
  if (x <= 0x3FFFFFFFFFFFFFFFL) {n = n + 2; x = x << 2;}
  if (x <= 0x7FFFFFFFFFFFFFFFL) {n = n + 1;}
  return n;
}

Дальше подсчитать несложно.


Ответ 2

Зачем терять точность? Ниже точная оценка числа единичных битов.
scanf("%u", &x);

x = (x & 0x55555555) + ((x >> 1) & 0x55555555);
x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
x = (x & 0x0F0F0F0F) + ((x >> 4) & 0x0F0F0F0F);
x = (x & 0x00FF00FF) + ((x >> 8) & 0x00FF00FF);
x = (x & 0x0000FFFF) + ((x >>16) & 0x0000FFFF);

printf("%d\n", x);

или
scanf("%u", &x);

x = x - ((x >> 1) & 0x55555555); 
x = (x & 0x33333333) + ((x >> 2) & 0x33333333); 
x = (x + (x >> 4)) & 0x0F0F0F0F; 
x = x + (x >> 8); 
x = x + (x >> 16); 
x &= 0x0000003F;

printf("%d\n", x);

А, я, кажется, неверно понял вопрос, имелось в виду сколько разрядов, считая от левой
единицы? Вот код, дающий количество ведущих нулей, 32 - n даст количество значащих
разрядов.
unsigned int x;
unsigned int n = 1;

scanf("%u", &x);

if (x == 0) return(32);

if ((x >> 16) == 0) {n = n + 16; x = x <<16;}
if ((x >> 24) == 0) {n = n + 8; x = x << 8;}
if ((x >> 28) == 0) {n = n + 4; x = x << 4;}
if ((x >> 30) == 0) {n = n + 2; x = x << 2;}
n = n - (x >> 31);

printf("%d\n", n);


Ответ 3

Можно сделать таблицу степеней 2 и проводить по ней бинарный поиск:
#include 

unsigned int table[32] = {
    0x00000001, 0x00000002, 0x00000004, 0x00000008,
    0x00000010, 0x00000020, 0x00000040, 0x00000080,
    0x00000100, 0x00000200, 0x00000400, 0x00000800,
    0x00001000, 0x00002000, 0x00004000, 0x00008000,
    0x00010000, 0x00020000, 0x00040000, 0x00080000,
    0x00100000, 0x00200000, 0x00400000, 0x00800000,
    0x01000000, 0x02000000, 0x04000000, 0x08000000,
    0x10000000, 0x20000000, 0x40000000, 0x80000000
};

int bits(unsigned int n) {
    int i=0, j=32, k=16;
    while(i>1;
    }
    return j;
}

int main() {
    int i;
    for(i=0;i<=0x100;i++)
        printf("0x%08X %d\n",i,bits(i));
    return 0;
}

P.S. А можно сделать бинарным поиском, но без таблицы :
int bits(unsigned int n) {
    int i=0, j=32, k=16;
    while(i>1;
    }
    return j;
}


Ответ 4

У многих процессоров есть подобная инструкция. Но писать под каждый процессор свою
реализацию не очень удобно.
В gcc есть встроенные функции __builtin_clz..., которые в случае поддержки процессора,
компилируются в соответствующие инструкции, иначе gcc сам генерит какой-то код.
В MSVS тоже есть нечто подобное, __lzcnt..., но в документации сказано, что разработчик
должен проверить поддержку этой операции процессором, иначе результат вызова непредсказуем.


Ответ 5

Исчо один вариант. С длинным сдвигом. Результат от 0 (для аргумента 0) до 32.
static int m[] = {0xffff0000, 0xff00, 0xf0, 0xC, 2};
int bits(unsigned long a) {
    int n=0;
    int k=16;
    int *m1 = m;
    while(k) {
        if(a&*m1++) {
            n += k;
            a >>= k;
        }
        k >>= 1;
    }
    return a?n+1:n; // (0..31) => (1..32) для a != 0
}

Любители сишных трюков могут перенести сдвиг k в заголовок while для большей кучерявости
(тогда нач. значение k будет 32)


Ответ 6

Странно, что никто не предложил использовать массив в качестве хэша:
int f(unsigned int n)
{
  static const int bitsByNum[16] = {1,1,2,2,3,3,3,3,4,4,4,4,4,4,4,4};

  int removed = 0;
  if (n > (unsigned int)0xFFFF)
  {
    removed += 16;
    n >>= 16;
  }
  if (n > (unsigned int)0xFF)
  {
    removed += 8;
    n >>= 8;
  }
  if (n > (unsigned int)0xF)
  {
    removed += 4;
    n >>= 4;
  }
  return removed + bitsByNum[n];
}


Ответ 7

Можно еще так попробовать, правда сомневаюсь, что это даст выиграш производительности,
разве что на числах очень большой разрядности: 
Предположим у вас есть N-битное число Y.
делаем xor между старшими N/2 битами числа и N/2 младшими (повторяем необходимое
количество раз)
результат каждой итерации имеет точность 0.5-1.0 от реального значения (т.е точность
2х будет 0.25 - 1, 0,0625 - 1).
Можно подобрать необходимое количество итераций, а что более важно можно будет прогнать
на реальных данных даную функцию и сравнивать ее точность с точностью исходного алгоритма
(приведенного в теле вопроса)


Ответ 8

А если так (на Java):
int d=Integer.toBinaryString(N).length();

Как называется по-русски распределение, плотность которого 1/x?

#математика #терминология #численные_методы #теория_вероятностей

Ричард Хэмминг опубликовал в 1970 году статью "On the Distribution of Numbers". В
этой статье он решил вопрос о том, как выглядит распределение мантисс всех представимых
чисел в арифметике с плавающей точкой по основанию b (мантисса, это число в интервале
[1/b;1] )

Решение выглядит очень просто - в пределе это распределение стремится к распределению
с плотностью

r(x) = 1 / ( x ln b )


По-английски с лёгкой руки Хэмминга это распределение называется Reciprocal Distribution.
Однако в учебниках на русском языке мне не удалось найти ничего об этом вероятностном
распределении.

Математические словари предлагают перевод 

reciprocal distribution => квантиль распределения


Этот перевод, как мне кажется, совсем не о том.

Как принято называть reciprocal distribution по-русски?
    

Ответы

Ответ 1

Инверсное распределение.

Сошлюсь на статью в Википедия

Метод Рунге-Кутты 4 порядка на java

#java #методы #численные_методы

Я только начала изучать java, и многие элементарные вещи для меня являются непонятными.

Прощу вас подсказать мне, как решать систему взаимозависимых уравнений методом Рунге-Кутты!

Полное задание:

Для решения полученной задачи Коши для системы первого порядка вида:

y'= f(t,y), y(0)=y0

использовать метод Рунге_Кутты 4-го порядка точности:

k1 = f(tn , yn)
 k2 = f(tn + h/2 , yn + hk1 / 2)
 k3 = f(tn + h/2 , yn + hk2 / 2)
 k4 = f(tn + h , yn + h*k3)

yn+1=yn+h*(k1 + 2*k2 + 2*k3 + k4) / 6    

На отрезке [0,5] с точным решением 

y1=cos(x)/(1+e2x)1/2,

y2=sin(x)/(1+e2x)1/2.

Для проверки правильности работы программы решить тестовую задачу из двух уравнений:

y1' = –y2 + y1(y12 + y22 – 1),

y2' = y1 + y2(y12 + y22 – 1),

На отрезке [0,5] с точным решением 

y1=cos(x)/(1+e2x)1/2,

y2=sin(x)/(1+e2x)1/2.

Я пытаюсь реализовать только метод!
Сложность заключается в том, что производные у1 и у2 зависят друг от друга.
И правильно ли  я сделала, что производную в точке вычисляю по у1 и у2

public static double y1=Math.cos(x)/Math.pow(1 + Math.pow(Math.E, 2 * x),0.5);
public static double y2=Math.sin(x)/Math.pow(1+Math.pow(Math.E,2*x),0.5);

//   dy1/dx
public static double derviY1(double x,double y10,double y20){

    return -y20+y10*(Math.pow(y10,2)+Math.pow(y20,2)-1);
}

//   dy2/dx
public static  double derviY2(double x ,double y1,double y2){
    return y1 + y2*(Math.pow(y1,2) + Math.pow(y2,2) - 1);
}


И почему-то вычисляя y20[i+1] и y10[i+1] они у меня остаются неизменными, хотя я
в программе изменяю данные, которые в них входят.

for (int i = 0; i < n - 1; i++) {
    x = i * h;

    k1 = h * derviY1(x, y10[i], y20[i]);
    m1 = h * derviY2(x, y10[i], y20[i]);

    k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2);
    m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);

    k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2);
    m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);

    k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3);
    m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);

    y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
    y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;

    System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " + y20[i]
+ " " + "|");


Полный код программы:

import java.*;
import java.lang.Math.*;
import static java.lang.System.out;

public class Test {

    static double x;
    //start
    private static int a = 0;
    // stop
    private static int b = 5;

    public static double y1 = Math.cos(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5);
    public static double y2 = Math.sin(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5);

    //   dy1/dx
    public static double derviY1(double x, double y10, double y20) {

        return -y20 + y10 * (Math.pow(y10, 2) + Math.pow(y20, 2) - 1);
    }

    //   dy2/dx
    public static double derviY2(double x, double y1, double y2) {
        return y1 + y2 * (Math.pow(y1, 2) + Math.pow(y2, 2) - 1);
    }


    public static void main(String[] args) {
        int n = 10;
        double h = (b - a) / n;
        double k1, k2, k3, k4, m1, m2, m3, m4;
        double[] y10 = new double[n]; //array of values y1
        double[] y20 = new double[n]; //array of values y2
        y10[0] = 1 / Math.sqrt(2);
        y20[0] = 0;

        // Computation by 4th order Runge-Kutta
        //update x
        for (int i = 0; i < n - 1; i++) {
            x = i * h;

            k1 = h * derviY1(x, y10[i], y20[i]);
            m1 = h * derviY2(x, y10[i], y20[i]);

            k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2);
            m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);

            k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2);
            m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);

            k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3);
            m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);

            y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;

            System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " +
y20[i] + " " + "|");
        }
    }
}

    

Ответы

Ответ 1

У вас h = 0. Почему? Ведь вы ясно написали, что h = (5-0)/10. А проблема вот в чём.
Ваше выражение (5-0)/10 сначала преобразуется в int, а потом уже в double. То есть
происходит всё как-то так (int)((5-0)/10). После этой операции мантисса отбрасывается
и получается полноценный 0.
Достаточно привести к double одну из переменных, как всё выражение будет обработано,
как для переменных с плавающей точкой. К примеру так:

double h = ((double)b - a) / n;


В этих преобразованиях есть более интересные моменты, о которых хотелось бы рассказать,
то это уже другая история.

Численное решение задачи теплопроводности

#python #численные_методы

Необходимо решить задачу теплопроводности на отрезке


При решении использовала явную схему

N = 10 # максимальное число шагов по х
K = 10 # максимальное число шагов по t
l = 1 # значение х на правой границе
h = l / N # шаг сетки по х
T = 1 # максимальное значение времени t на правой границе
t = T / K # шаг сетки по времени

# зададим сетку 
x_i = np.arange(0, N, h) # значения в узлах по х
t_j = np.arange(0, K, t) # значение в узлах по t
r_j = len(t_j) # количество узлов по t
r_i = len(x_i) # количество узлов по x
w_h_t = np.zeros([r_i, r_j]) # итоговая сетка размером x_i*t_j

# зададим значение функции входящей в начальное уравнение
x = 0
def f(x):
    return np.sin(x)

# граничные условия
ux_0 = 1 # граничное условие на левом конце при x=0
ut_0 = np.cos(x_i) # граничное условие при t=0

# найдем значения на нулевом слое при t=0 ut_0 = np.cos(x_i)
w_h_t[0] = np.cos(x_i)

# найдем значения w_h_t на первом и последующих слоях
const = t / (h**2) 
for j in range(1, len(x_i)-1):
    for i in range(len(w_h_t[j])-1):        
        w_h_t[j+1, i] = w_h_t[j, i] + const * (w_h_t[j,i] - 2*w_h_t[j,i] + w_h_t[j,
i-1]) + t*f(x_i[j])
        w_h_t[j+1, 0] = 1
        w_h_t[j+1, len(w_h_t[i])-1] = w_h_t[j+1, len(w_h_t[i])-2] + h * t_j[j+1]


plot_ = np.arange(0,len(w_h_t)-1,1)
for y in plot_:
    plt.plot(x_i, w_h_t[y])


В результате получила вот такой график зависимости Х от рассчитанного значения функции
в узлах сетки 



Мне не понятно насколько неверно мое решение. Возникли проблемы с поиском частного
решения, wolfram выдал u(x) = 1.54x+1+sinx, что мне кажется не верным, а самой решить
не получилось. В учебнике Филиппова похожих примеров не нашлось, в сети ничего достаточно
подробного, чтобы разобраться в решении не нашла. Подскажите где можно найти как решать
аналитически такое уравнение и насколько неверно мое решение? В чем ошибка? И как вообще
проверяют на корректность решения таких задач, кроме как сравнения с аналитическим
решением?



В общем и целом разобралась. Решила дополнить свой вопрос отредактированным решением,
возможно, кому-то пригодится.
Численное решение правильно найти так и не удалось, но находится оно методом Фурье(разделение
переменных). 
Итоговый график:
 

Рабочий код:

N = 10 # максимальное число шагов по х
K = 500 # максимальное число шагов по t
l = 1 # значение х на правой границе
h = l / N # шаг сетки по х
T = 1 # максимальное значение времени t на правой границе
t = T / K # шаг сетки по времени


# зададим сетку 
x_i = np.arange(0, 1, h) # значения в узлах по х
t_j = np.arange(0, 1, t) # значение в узлах по t
r_j = len(t_j) # количество узлов по t
r_i = len(x_i) # количество узлов по x
w_h_t = np.zeros([r_j, r_i]) # итоговая сетка размером x_i*t_j


# зададим значение функции входящей в начальное уравнение
x = 0
def f(x):
    return np.sin(x)

# граничные условия
ux_0 = 1 # граничное условие на левом конце при x=0
ut_0 = np.cos(x_i) # граничное условие при t=0

# найдем значения на нулевом слое при t=0 ut_0 = np.cos(x_i)
w_h_t[0] = np.cos(x_i)

# найдем значения w_h_t на первом и последующих слоях
const = t / (h**2) 
for j in range(len(w_h_t) - 1):
    for i in range(len(w_h_t[j]) - 1):        
        w_h_t[j + 1, i] = w_h_t[j, i] + const* (w_h_t[j, i+1] - 2 * w_h_t[j, i] +
w_h_t[j, i - 1]) + t*f(x_i[i])
        w_h_t[j + 1, 0] = 1
        w_h_t[j + 1, len(w_h_t[i])-1] = w_h_t[j + 1, len(w_h_t[i])-1] + h

    

Ответы

Ответ 1

К сожалению не специалист именно в теплопроводности, но могу отметить несколько точно
проблемных моментов:


В настоящий момент шаг по сетке времени не связан с шагом по сетке в пространстве.
В явных схемах - это ведет к нестабильности. Нужно соблюдать критерий CFL. В вашем
случае, если я правильно помню термодинамику, Δt < CFL * χ * (Δx)². То есть Δt < (Δx)²
(χ - radiative diffusion в вашем уравнении = 1). В вашей программе это явно не соблюдается
ни по конкретным значениям, ни по алгоритму. Решение: привяжите K и N друг к другу.
По вашему графику очень сложно понять что происходит. Если начертить каждый шаг времени
на отдельном графике - то видно, что численное решение расходится. В частности, это
видно и на вашем чертеже (e121). Объясняется как минимум пунктом 1, и, возможно, еще
какими-то проблемами в программе. Кроме того, решения на промежуточных шагах выглядят
ну очень подозрительно.
Посмотрите значения x_i и t_j (print(x_i)) - сейчас это массивы от 0 до 9.9  с шагом
0.1. Судя по условию, перепутан шаг, количество шагов и правая граница.


Теперь по проверке решений. Формально, такой код можно проверить через Method of
Manufactured Solutions, но для вашего задания это, конечно, слишком. 

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


выполняются ли начальные и конечные условия
анализ сходимости (убавили шаг по времени - стало лучше? убавили шаг дискретизации
в пространстве - стало лучше?)
не появляется ли энергия из ниоткуда (для замкнутых систем)?
совпадает ли решение с решением из коммерческого\open-source симулятора?
если решить уравнение другим численным методом (например, методом конечных элементов),
получим ли мы примерно тоже самое?

МНК алгоритм

#алгоритм #cpp #численные_методы

Метод наименьших квадратов
Задача написать алгоритм для нахождения коэффициентов полинома N степени с помощью
МНК по M точек. То есть даны 10 точек, и сказано аппроксимировать для полинома 3-й
степени: 
y=ax^3 + a1x^2 + a2x+a3

Для этого полинома надо найти a, a1, a2, a3; если для 1 и 2 степени можно написать
вручную, то даже для 3-й степени придётся решать систему 4-х уравнений с 4 неизвестными.
Вот тут можно посмотреть, как выглядит аппроксимация для разных функций с помощью мнк.
Как запрограммировать алгоритм, который будет составлять СЛОУ и решать её?
Язык c++.
Обновление
Как в составлении так и в решении. Начал решать пока что такая система получается
для полинома N степени. Ex сумма от 1 до M (кол-во точек).
a1*Ex^n     + a2*Ex^n - 1  + ... + an*Ex^0=Eyx^0

a1*Ex^n + 1 + a2*Ex^n      + ... + an*Ex^1=Eyx^1


a1*Ex^2n    + a2*Ex^2n - 1 + ... + an*Ex^n=Eyx^n

Даже после того как перепишу все известные Ex,Ey,Exy... в массив, то я незнаю как
решать такую СЛАУ
Решение
Прошу меня простить за кривые отступы.
Функция которая находит коэф. полинома N степени
    vector* mnk(vector> pts,int n,int n2)  
    {
        //матрица хранящая неизвестные при нужных нам коэффициентах 
        vector> linSys(n + 1);
double tmp=0;
vector subsid(3*n + 2);//почему 3n + 2? да потому что http://mathhelpplanet.com/static.php?p=onlayn
- mnk - i - regressionniy - analiz
for (int i=0;i<2*n + 1;i ++ )//этот цикл формирует первую половину (точнее 2/3) списка
величин x^0 x^1...x^2n + 1
{
    for (int j=0;j*b=new vector;
b - >resize(n + 1);
for(int i = 0; i < n  +  1; i ++ )
{
    linSys[i].resize(n + 1);
    (*b)[i]=subsid[2*n + 1 + i];
    for(int j=0;j*x=new vector;
x - >resize(n + 1);
//это решает СЛАУ и пишет результ в x
Gauss(linSys,*b,x,n + 1);
return x;

}
Функция которая решает СЛАУ 
#include 
#include 
#include 
#include 
//результат записывает в x, СЛАУ представима как A*x=B
void Gauss(vector> A, vector b, vector *x, int n)
    {
std::stack > swaps;
int i, j, k, t;
double kof, s;
for (i = n  -  1; i > 0;  -- i) {
    for (t=i, j=i - 1; j >= 0;  -- j) {//Поиск максимума в столбце
        if (fabs(A[i][t]) < fabs(A[i][j])) {
            t = j;
        }
    }
    if (A[i][t] == 0.0) {
        return;
    }
    if (t != i) {
        for (k = n  -  1; k >= 0;  -- k) {
            std::swap(A[k][t],A[k][i]);
        }
        swaps.push(std::pair(i,t));
    }
    for (j=i - 1; j>=0;  -- j) {//
        kof = A[j][i] / A[i][i];//
        A[j][i] = 0.0;//
        b[j]  - = b[i] * kof;
        for (k = i  -  1; k >= 0;  -- k) {
            A[j][k]  - = A[i][k] * kof;
        }
    }
}
for (i = 0; i < n;  ++ i)
{
    s = 0.0;
    for (j = i  -  1; j >= 0;  -- j) {
        s  + = A[i][j] * (*x)[j];
    }
    (*x)[i] = (b[i]  -  s) / A[i][i];
}
while (!swaps.empty()) {
    //int j=swaps.top().first;
//  int y=swaps.top().second;
    std::swap((*x)[swaps.top().first], (*x)[swaps.top().second]);
    swaps.pop();
}

}    

Ответы

Ответ 1

Задача о полиномиальной регрессии (аппроксимации табличной функции) требует большой
аккуратности, поскольку при обычных методах решения расчёт главного определителя СЛАУ
(определителя Вандермонда) приводит к вычитанию близких чисел и в итоге снижает относительную
погрешность вычислений, что проявляется уже для полиномов 8-10 порядка. 

От этих недостатков свободны методы разложения в ряды по ортогональным функциям (например,
разложение в ряды Фурье по косинусам), и хотелось бы так же просто обращаться с полиномами.
Однако классические ортогональные полиномы (Лежандра, Чебышёва и пр.) в дискретном
варианте ортогональны только на специально подобранных неравномерных сетках.

Удачный выход из такой ситуации даёт малоизвестный метод ортогональной полиномиальной
регрессии (Orthogonal Polynomial Curve Fitting, Jeff Reid), в котором сначала по абсциссам
заданных точек конструируются ортогональные полиномы, а затем с помощью МНК рассчитываются
коэффициенты при этих полиномах. 

1. ТЕОРИЯ

1.1. Постановка задачи

Даны n точек { (xi,yi), i = 0, 1 ... n-1 }.
Найти полином g(x) = a0 + a1x + ... + amxm с минимальной суммой ∑i = 0, 1 ... n-1
(yi - g(xi))2 среди всеx возможных полиномов n-го порядка.  

Решение будем искать в виде g(x) = b0 p0(x) + b1 p1(x) + ... + bm pm(x),
где pj(x), j = 0, 1, ... m - семейство ортогональных полиномов, для которых должны
выполняться рекуррентные соотношения

pj+1(x) = (x-Aj+1) pj(x) -  Bj pj-1(x), p0(x) = 1, p-1(x) = 0, (1)

и условия ортогональности

∑i = 0, 1 ... n-1 pj(xi) pk(xi) = 0  при  j ≠ k. (2)

1.2. Построение ортогональных полиномов.

Для j = 0 формулы (1) и (2) дают:  p1(x) = x - A1,  ∑ (xi - A1) = 0, откуда

A1 = ∑xi / n.

Если для j = 0 ... k ортогональные полиномы pj уже известны, то следующий полином
pk+1 должен быть ортогонален каждому из них:  ∑ pj(xi) pk+1(xi) = 0.  Используя рекурррентные
соотношения (1), получаем: 

∑ xi pj (xi) pk (xi) - Ak+1∑ pj (xi) pk (xi) - Bk∑ pj(xi) pk-1 (xi) = 0,  (3) 

xpj(x) = pj+1(x) + Aj+1pj(x) +  Bj pj-1(x),  (1')

Из условий (2) и (1') следует, что для j = 0 ... k-2 все три слагаемых в соотношениях
(3) нулевые. При j = k - 1 условия ортогональности уничтожают второе слагаемое в формуле
(3) и два младших слагаемых в (1'), поэтому

Bk = ∑ p2k(xi) / ∑ p2k-1(xi).  (4) 

В случае j = k условия ортогональности обнуляют третье слагаемое в формуле (3), поэтому

Ak+1 = ∑ xip2k(xi) / ∑ p2k(xi).  (5)  

1.3. Построение полинома регрессии

Метод наименьших квадратов для полинома
g(x) = b0p0(x) + b1p1(x) + ... + bmpm(x)
предполагает вычисление коэффициентов bj, j = 0 ... m, по методу наименьших квадратов,
т.е. минимизацию функции невязки
F(b0, b1 ... bm) = ∑ (yi - b0p0(xi) - b1p1(xi) - ... - bmpm(xi))2
по этим коэффициентам, рассматриваемым в качестве переменных. При этом в точке минимума
должны выполняться необходимые условия экстремума, т.е. равенство нулю частных производных
F'bk(b0, b1 ... bm):

-2∑ (yi - b0p0(xi) - b1p1(xi) - ... - bmpm(xi)) pk = 0,  k = 0 ... m.

Разбивая каждое из этих уравнений на слагаемые и используя условия ортогональности
(2), получим:
∑ yipk(xi) - bk ∑ pk2(xi) = 0,  k = 0 ... m,

bk = ∑ yipk(xi) / ∑ pk2(xi),  k = 0 ... m. (6)

Приводя полином g(x) к виду
g(x) = a0 + a1x + ... + amxm,
получим:
aj = ∑k=j...m bk pkj,  j = 0 ... m. (7)

Формулы (6-7) позволяют построить алгоритм расчёта коэффициентов 
aj  полинома g(x), если ортогональные полиномы известны.

2. АЛГОРИТМ

2.1. Значения ортогональных полиномов в узлах сетки

Значения полиномов Pji = pi(xi) в узлах сетки можно вычислить по формулам (1,4-5):

P0i = 1, 
P1i = xi -  Q0 / S0, 
Pj+1,i = (xi -  Qj / Sj) Pj,i - (Sj / Sj-1) Pj-1,i для j > 0,  

где Sj = ∑i Pj,i2, Qj = ∑i xi Pj,i2.

Для полинома третьего порядка (m=3) алгоритм имеет вид:

P0i = 1;  S0 = m+1;  Q0 = ∑ xi,

P1i = xi - Q0 / S0;  S1 = ∑ (P1i)2;  Q1 = ∑ xi (P1i)2,

P2i = (xi - Q1 / S1) P1i - S1 / S0;  S2 = ∑ (P2i)2;  Q2 = ∑ xi (P2i)2,

P3i = (xi - Q2 / S2) P2i - S2i / S1i;  S3 = ∑ (P3i)2;  Q3 = ∑ xi (P3i)2.

2.2. Коэффициенты ортогональных полиномов

Рекуррентные формулы (1,4,5) также позволяют вычислить все коэффициенты сkj ортогональных
полиномов вида pk(x) = сk0 + сk1 x + ... + сkj xj +... + сkm xm, которые в общем случае
представляют собой алгебраическую сумму трёх величин:

сk+1,j = сk,j-1 - (Qk / Sk) сk,j - (Sk / Sk-1) сk-1,j для k = 0, 1 ... m-1, j = 0,
1, ... , k, 

сk,k = 1 для  k = 0, 1 ... m.

Первое слагаемое не следует учитывать при j = 0, третье - при j = k.

Для полинома третьего порядка (m=3):

с00 = 1,  с10 = - Q0 / S0;  с11 = 1,

с20 = - (Q2 / S2) с10 - (S1 / S0) с00;  
с21 = с10 - (Q2 / S2) с11 - (S1 / S0) с01; 
с22 = 1,

с30 = - (Q3 / S3) с20 - (S2 / S1) с10; 
с31 = с20 - (Q3 / S3) с21 - (S2 / S1) с11; 
с32 = с21 - (Q3 / S3) с22 - (S2 / S1) с12;  
с33 = 1.

2.3. Коэффициенты полинома регрессии.

Коэффициенты bk в разложении g(x) = b0 p0(x) + b1 p1(x) + ... + bm pm(x) вычисляются
в соответствии с (6) по формулам

bk = ∑i yi Pki / Sk,  k = 0 ... m.  

Для полинома третьего порядка:

b0 = ∑i yi P0i / S0, 
b1 = ∑i yi P1i / S1, 
b2 = ∑i yi P2i / S2, 
b3 = ∑i yi P3i / S3. 

Коэффициенты aj в разложении g(x) = a0 + a1 x + ... + am xm вычисляются по формуле
aj = ∑ k=j...m bk сkj,  j = 0 ... m. (7)

Для полинома третьего порядка (m=3):

a0 = b0 с00 + b1 с10 + b2 с20 + b3 с30,
a1 = b1 с11 + b2 с21 + b3 с31,
a2 = b2 с22 + b3 с32,
a3 = b3.

2.4. Об алгоритме

Разработанный алгоритм в общем виде применим для аппроксимации заданного массива
точек полиномом любого порядка при условии достаточности исходных данных и вычислительных
средств.
На примере полинома третьего порядка показано, что никакие СЛАУ при этом решать не
требуется.

3. ТЕСТИРОВАНИЕ АЛГОРИТМА

Тестирование алгоритма проведено с использованием пакета MathCAD.
Для тестирования были выбраны 10-точечные массивы x и y и полиномиальная модель третьего
порядка (см. скриншот).



Тестирование подтвердило правильность алгоритма.

Сравнение результатов расчёта с функцией regress() пакета MathCAD показало, что при
одинаковой разрядности расчётов (определяемой исходными данными) алгоритм ортогональной
полиномиальной регрессии точнее.  

4. ПРОГРАММА НА PHP

Программа сoдержит класс Ortho, конструктор которого принимает исходные массивы $x,
$y и порядок модели $m и рассчитывает регрессионный полином.
Для выдачи рассчитанных величин предусмотрены следующие методы:
getValues() - для массива значений ортогональных полиномов в точках x;
getSums() - для сумм, рассчитанных по этому массиву;
getAllCoeffs() - для массивов коэффициентов (с - массивы коэффициентов ортогональных
регрессионных полиномов; b - массив для расчёта регрессии по значениям ортогональных
полиномов; a - коэффициенты при степенях регрессионного полинома);
getRegress() - для значений регрессионного полинома и их невязок по отношению к массиву $y.

/**
*   Polinomial regression via orthogonal regression polynomials.
*   Author: Yury Negometyanov 
*   21.03.2017 
*/

class Ortho
{
    private $x;
    private $y;
    private $m;
    private $values;
    private $sums;
    private $coeffs;

    public function __construct($xx, $yy, $mm)
    {
        $this->x = $xx;
        $this->y = $yy;
        $this->m = $mm;
        $this->calcValuesSums();
        $this->calcAllCoeffs();
    }

    public function subtract($ar, $subtr)
    {
        if(gettype($subtr) == "array") {
            return array_map(function($val1, $val2){ return $val1 - $val2;}, $ar,
$subtr);      
        } else { 
            return array_map(function($val)use($subtr){ return $val - $subtr;}, $ar);
      
        }
    }

    public function multi($ar, $mult)
    {
        if($mult == "square") {
            return array_map(function($a){return $a * $a;}, $ar);
        } elseif(gettype($mult) == "array") {
            return array_map(function($a, $b){return $a * $b;}, $ar, $mult);
        } else {
            return array_map(function($val)use($mult){ return $val * $mult;}, $ar);     
        }
    }

    private function calcValuesSums()
    {
        for ($j = -1; $j < $this->m ; $j ++) {
            if(!isset($v)) {
                $v = [array_fill(0, count($this->x), 1.0)];
                $s = [(float)count($this->x)];
                $q = [array_sum($this->x)];
            } else {
                $v[] = $this->subtract($this->x, end($q)/end($s)); 
                if($j > 0)  {
                    $v[$j+1] = $this->subtract( $this->multi(end($v), prev($v)), 
                                                $this->multi(prev($v), end($s)/prev($s))); 
                }
                $sq = $this->multi(end($v), "square");
                $s[] = array_sum($sq);
                $q[] = array_sum($this->multi($this->x, $sq));
            }
        }
        $this->values = $v;
        $this->sums = ['∑P²'=>$s, '∑xP²'=>$q];
    }

    private function calcOrthoCoeffs()
    {
        $prev = null;
        $ab = array_map(function($ss, $qq)use(&$prev)
        {
            if(is_null($prev)) $bb = 0;
            else $bb = $ss / $prev;
            $prev = $ss;
            return [$qq / $ss, $bb];
        }, reset($this->sums), next($this->sums));
        foreach ($ab as $k => list($aa, $bb)) {
            if(!isset($с)) $с = [[1.0]];
            $с[] = array_fill(0, $k+2, 1.0);
            for ($j=0; $j <= $k; $j++) { 
                $с[$k+1][$j] = ($j==0 ? 0.0 : $с[$k][$j-1]) 
                                - $с[$k][$j] * $aa 
                                - (($j==$k) ? 0.0 : $с[$k-1][$j] * $bb);
            }
        }
        return $с;
    }

    private function calcAllCoeffs()
    {
        $s = reset($this->sums);
        $b = [];
        foreach ($this->values as $key => $v) {
            $b[] = array_sum($this->multi($this->y, $v))/$s[$key];
        }
        $c = $this->calcOrthoCoeffs();
        $m = count($b) - 1;
        $a = array_fill(0, $m+1, 0.0);
        foreach ($a as $j => &$aj) {
            for ($k = $j; $k <= $m; $k++) {
                $aj += $b[$k] * $c[$k][$j]; 
            }
        }
        $this->coeffs = ['c' => $c, 'b'=>$b, 'a'=>$a];
    }

    private function calcRegress($opt ='a')
    {
        switch ($opt) {
            case 'a':
                $rev = array_reverse($this->coeffs['a']);
                $p = [];
                foreach($this->x as $xx) {
                    $sum = 0.0;
                    foreach ($rev as $value) {
                        $sum *= $xx;
                        $sum += $value;
                    }
                    $p[] = $sum;
                }
                break;

            case 'b':   
                foreach ($this->values as $key => $arr) {
                    if(!isset($p)) $p = $this->multi($arr, $this->coeffs['b'][$key]);
                    else $p = $this->subtract($p, $this->multi($arr, - $this->coeffs['b'][$key]));
                }
                break;

            default:
                $p = [];
                break;
        }
        return $p;
    }

    public function getValues() 
    {
        return $this->values;
    }

    public function getSums() 
    {
        return $this->sums;
    }

    public function getAllCoeffs() 
    {
        return $this->coeffs;
    }

    public function getRegress() 
    {
        $values_a = $this->calcRegress('a');
        $values_b = $this->calcRegress('b');
        return [
            'regress_a'=>$values_a,
            'regress_b'=>$values_b,
            'discrep_a'=>array_sum($this->multi($this->subtract($this->y, $values_a),
"square")), 
            'discrep_b'=>array_sum($this->multi($this->subtract($this->y, $values_b),
"square"))
        ]; 
    }
}

function print_m($text, $arr, $level=0){
    $space = str_repeat(" ", $level++);
    echo "$space$text";
    $flag = false;
    foreach($arr as $value) $flag = $flag || (gettype($value)=="array");
    foreach($arr as $key => $value) {
        if(gettype($value) != "array") {
            echo $flag ? "
$key => $value" : " $value";
        } else {
            print_m("
$space$key => [", $value, $level);
            echo " ]";
        }
    }
    $level--;
}   

$x = [0.00000,  3.36588,  3.63719, 0.56448,  -3.02721,  -3.83570, -1.11766,  2.62795,
  3.95743,  1.64847];
$y = [3.95610, 74.84479, 89.44289, 6.46668, -14.53888, -34.55881,  1.70531, 43.80101,
109.12940, 18.81613];

/*
$x = [];
for ($i=0; $i < 1000; $i++) { 
    $x[] = (float)5*sin($i);
}

$y = [];    
foreach ($x as $xx) {
    $y[] = (float)sin($xx/2+1);
}

$n = count($x);
$ort = new Ortho($x, $y, $m = 80);
$discrep = $ort->getRegress();

echo "Sample size = $n  Model Order = $m";
echo "
Discrepancy via Powers is: {$discrep['discrep_a']}";
echo "
Discrepancy via Ortogonal Regression Polynomials is: {$discrep['discrep_b']}";
*/

$ort = new Ortho($x, $y, $m=3);
print_m("Issue Data:", ['Array_x'=>$x, 'Array_y'=>$y, 'Model Order'=>$m]);
$values = $ort->getValues();
print_m("

The Orthogonal Regression Polynomials' Values:", $values);
$sums = $ort->getSums();
print_m("

Calculated Sums:", $sums);
$coeffs = $ort->getAllCoeffs();
print_m("

Calculation of Polynomials' Coefficients:", $coeffs);
$discrep = $ort->getRegress();
print_m("

Regression via Powers (a) & via Orthogonal Regression Polynomials
(b) :", $discrep);


5. РЕЗУЛЬТАТЫ.

Результаты на тестовой выборке:


Issue Data: 
Array_x => [ 0 3.36588 3.63719 0.56448 -3.02721 -3.8357 -1.11766 2.62795 3.95743 1.64847 ] 
Array_y => [ 3.9561 74.84479 89.44289 6.46668 -14.53888 -34.55881 1.70531 43.80101 109.1294 18.81613 ]
Model Order => 3

The Orthogonal Regression Polynomials' Values: 
0 => [ 1 1 1 1 1 1 1 1 1 1 ] 
1 => [ -0.782083 2.583797 2.855107 -0.217603 -3.809293 -4.617783 -1.899743 1.845867 3.175347 0.866387 ] 
2 => [ -7.2281734230875 2.8073623836192 4.6030924341892 -7.1264829728247 3.0992779145833 8.9585998726813 -5.5494580486592 -1.3320551385999 6.9121153510662 -5.1442783729677 ] 
3 => [ 3.6796621840027 -2.8171823354639 3.1956859234941 -3.0255967835539 8.7399448097287 -12.367028020194 15.20316889478 -12.281110609376 12.297473127759 -12.625017191178 ]

Calculated Sums: 
∑P² => [ 10 69.17098425001 328.7768235086 962.75401849569 ] 
∑xP² => [ 7.82083 -27.512890311149 -1.7129366948718 250.39738237352 ]

Calculation of Polynomials' Coefficients: 
c => [  
 0 => [ 1 ]  
 1 => [ -0.782083 1 ]  
 2 => [ -7.2281734230875 -0.38433110143359 1 ]  
 3 => [ 3.6796621840027 -11.983278954678 -0.37912107270775 1 ]  
 4 => [ 20.209167857339 7.9217601883617 -14.812965853531 -0.63920555697172 1 ] ] 
b => [ 29.906462 15.897655646369 2.3791287087584 1.000001267701 ] 
a => [ 3.9560877250835 2.9999883433859 2.0000071554385 1.000001267701 ]

Regression via Powers (a) & via Orthogonal Regression Polynomials (b) : 
regress_a => [ 3.9560877250835 74.844667502068 89.443009253453 6.4666635861522 -14.538829417751 -34.558843534167 1.7053151756382 43.801163134951 109.12933595098 18.816050623593 ] 
regress_b => [ 3.9560877250835 74.844667502068 89.443009253453 6.4666635861522 -14.538829417751 -34.558843534167 1.7053151756382 43.801163134951 109.12933595098 18.816050623593 ]
discrep_a => 6.7210313148693E-8
discrep_b => 6.7210313149619E-8

Результаты работы программы на тестовой выборке соответствуют ожидаемым.
На выборке в n=1000 точек при порядке модели 80 получено:


Sample size = 1000  Model Order = 80
Discrepancy via Powers is: 0.3261862533562
Discrepancy via Ortogonal Regression Polynomials is: 4.2188699484271E-26


Таким образом, расчёт регрессии через ортогональные полиномы более устойчив по отношению
к ошибкам округления.

6. ВЫВОДЫ.  


Алгоритм можно рекомендовать как удачную альтернативу традиционным алгоритмам для
полиномиальной регрессии произвольного порядка.
Исходные данные и полученные результаты могут быть использованы для отладки на других
языках программирования.

Алгоритм для гравитационной задачи N тел

#алгоритм #физика #численные_методы #моделирование

К своему стыду очень слаб в численных методах.
Суть задачи в том, чтобы смоделировать движение тел с течением времени под воздействием
силы гравитации друг на друга в количестве >= 3.

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

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

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


Например тела А и В имеют такое влияние друг на друга, что их атом времени составляет
5 секунд. При этом атом тел А и С составляет 1 / 1000000 секунды. Во время течения
атома A-B из-за гравитационного маневра тело А начало стремительное движение в сторону
тела В. В этот момент надо менять атом времени между A и В, иначе они пролетят мимо
друг друга. Получается нужна реактивность: при изменении стейта тела необходимо менять
атомы времени всех его связей с другими телами. Я правильно рассуждаю?
Должна ли при регулировке атома времени учитываться скорость сближения объектов?
Ведь два объекта могут быть достаточно далеко друг от друга чтобы не оказывать практически
никакое влияние друг на друга, но двигаться на сближение с огромной скоростью. Атом
времени между ними может оказаться настолько большим, что они сначала пролетят мимо
друг друга, и только потом высчитается их гравитационное влияние.


Сломал голову придумывая алгоритм реализации. Может быть есть у кого-нибудь пример
решения такой задачи: с радостью бы посмотрел код. А может-быть я услышу волшебное
"почитай про метод Вальфгауна-Штауца" и уйду гуглить в правильном направлении.

P.S. Для того чтобы немного поднять интерес прикреплю GUI существующего проекта.
Практической ценности к вопросу не несет никакого:)

https://makarov-andrey.github.io/n-body-problem/
    

Ответы

Ответ 1

Задача трех тел, есть частный случай задачи Коши для обыкновенных дифференциальных
уравнений: дан вектор начального состояния и функция который вычисляет производную
это вектора по времени (т.н. "функция правой части"). 
(Гуглить "численное решение задачи Коши".)

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

Типичный подход для оценки шага по времени (кроме привлечения физических соображений)
такой:


из начального состояния делаем шаг длинной t запоминаем
результат.  
из начального состояния делаем два шага длинной t/2
запоминаем результат. 
если результаты (1) и (2)  отличаются меньше чем заданный пользователем параметр
"точность на шаге", то принимаем этот шаг если больше то возвращаемся к начальному
состоянию и пытаемся сделать шаг меньшего размера. 
Насколько изменить шаг при пересчете шага, и как выбрать начально приближение для
следующего шага - зависит от выбранного метода решения. Для метода Ньютона, который
ты реализовал, чаще всего, делят или умножают длину шага на 2.


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

Важное замечание: задача трех - задача с разбеганием траекторий. Т.е. невозможно
заранее сказать с какой точностью нужно было посчитать решение для момента t1, чтобы
на его основании можно было получить решение для t2 с заданной точностью, пока не посчитаем
t2. Т.е. невозможно заранее сказать какая "точность на шаге" нужна для того чтобы получить
решение в момент t с заданной точностью. Типичный подход: задачу Коши решают несколько
раз, решают несколько раз, увеличивая "точность на шаге", пока результаты решений в
конечной точке по времени не будут отличаться друг от друга меньше чем заданная точность
интегрального решения.

Если цель изучения численных методов, рекомендую взять готовую библиотеку решения
задачи Коши. 
Для C++/C/fortran рекомендую:
 SUNDIALS https://computation.llnl.gov/projects/sundials 
Несколько самых простых методов, также реализованы в boost https://www.boost.org/doc/libs/1_64_0/libs/numeric/odeint/doc/html/index.html
(она попроще в использовании, но и методы послабее.)

Ключевые слова для поиска библиотек для других языков: метод Ньютона, метод Эйлера,
метод Розенброка, Stiffl, DASKR. 

Если цель изучить численные методы и построение солверов ОДУ, рекомендую книжку Э.Хайрер,
С.Нёрсет Г.Ваннер "Решение обыкновенных дифференциальных уравнений". (Поскольку система
уравнений в задаче трех тел - жесткая, тебе потребуются методы из начала 2-го тома:
"Жесткие и дифференциально-алгебраические задачи",  но без первого тома его читать
бесполезно.)

Интерактивные Графики в Python: как реализовать аналог функции Manipulate из Mathematica

Добрый день! Недавно начал изучать Python, решая задачу по численным методам. Хотел узнать, можно ли как-то в Python реализовать аналог функции Manipulate как в Mathematica, чтобы была возможность варьировать данные, чтобы это было видно на графике. Привел скрин ниже:

Буду рад любой помощи и подсказке. Спасибо!

Ответ

Можно использовать ipywidgets.interact (чтобы ползунки/график появились, ячейку возможно пару раз исполнить нужно)
import matplotlib.pyplot as plt import numpy as np from ipywidgets import interact
def f(a, k): x = np.arange(11) plt.plot(x, x**k - 1) plt.plot(x, a * (x - 1))
interact(f, a=(0,5), k=(0,6))

Можно записать этот код, используя синтаксис для аннотаций и декоратора:
@interact def f(a: (0,5), k: (0,6)): x = np.arange(11) plt.plot(x, x**k - 1) plt.plot(x, a * (x - 1))

Метод Рунге-Кутты 4 порядка на java

Я только начала изучать java, и многие элементарные вещи для меня являются непонятными.
Прощу вас подсказать мне, как решать систему взаимозависимых уравнений методом Рунге-Кутты!
Полное задание:
Для решения полученной задачи Коши для системы первого порядка вида:
y'= f(t,y), y(0)=y0
использовать метод Рунге_Кутты 4-го порядка точности:
k1 = f(tn , yn) k2 = f(tn + h/2 , yn + hk1 / 2) k3 = f(tn + h/2 , yn + hk2 / 2) k4 = f(tn + h , yn + h*k3)
yn+1=yn+h*(k1 + 2*k2 + 2*k3 + k4) / 6
На отрезке [0,5] с точным решением
y1=cos(x)/(1+e2x)1/2
y2=sin(x)/(1+e2x)1/2
Для проверки правильности работы программы решить тестовую задачу из двух уравнений:
y1' = –y2 + y1(y12 + y22 – 1),
y2' = y1 + y2(y12 + y22 – 1),
На отрезке [0,5] с точным решением
y1=cos(x)/(1+e2x)1/2
y2=sin(x)/(1+e2x)1/2
Я пытаюсь реализовать только метод! Сложность заключается в том, что производные у1 и у2 зависят друг от друга. И правильно ли я сделала, что производную в точке вычисляю по у1 и у2
public static double y1=Math.cos(x)/Math.pow(1 + Math.pow(Math.E, 2 * x),0.5); public static double y2=Math.sin(x)/Math.pow(1+Math.pow(Math.E,2*x),0.5);
// dy1/dx public static double derviY1(double x,double y10,double y20){
return -y20+y10*(Math.pow(y10,2)+Math.pow(y20,2)-1); }
// dy2/dx public static double derviY2(double x ,double y1,double y2){ return y1 + y2*(Math.pow(y1,2) + Math.pow(y2,2) - 1); }
И почему-то вычисляя y20[i+1] и y10[i+1] они у меня остаются неизменными, хотя я в программе изменяю данные, которые в них входят.
for (int i = 0; i < n - 1; i++) { x = i * h;
k1 = h * derviY1(x, y10[i], y20[i]); m1 = h * derviY2(x, y10[i], y20[i]);
k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2); m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);
k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2); m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);
k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3); m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);
y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6; y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;
System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " + y20[i] + " " + "|");
Полный код программы:
import java.*; import java.lang.Math.*; import static java.lang.System.out;
public class Test {
static double x; //start private static int a = 0; // stop private static int b = 5;
public static double y1 = Math.cos(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5); public static double y2 = Math.sin(x) / Math.pow(1 + Math.pow(Math.E, 2 * x), 0.5);
// dy1/dx public static double derviY1(double x, double y10, double y20) {
return -y20 + y10 * (Math.pow(y10, 2) + Math.pow(y20, 2) - 1); }
// dy2/dx public static double derviY2(double x, double y1, double y2) { return y1 + y2 * (Math.pow(y1, 2) + Math.pow(y2, 2) - 1); }
public static void main(String[] args) { int n = 10; double h = (b - a) / n; double k1, k2, k3, k4, m1, m2, m3, m4; double[] y10 = new double[n]; //array of values y1 double[] y20 = new double[n]; //array of values y2 y10[0] = 1 / Math.sqrt(2); y20[0] = 0;
// Computation by 4th order Runge-Kutta //update x for (int i = 0; i < n - 1; i++) { x = i * h;
k1 = h * derviY1(x, y10[i], y20[i]); m1 = h * derviY2(x, y10[i], y20[i]);
k2 = h * derviY1(x + h / 2, y10[i] + k1 / 2, y20[i] + k1 / 2); m2 = h * derviY2(x + h / 2, y10[i] + m1 / 2, y20[i] + m1 / 2);
k3 = h * derviY1(x + h / 2, y10[i] + k2 / 2, y20[i] + k2 / 2); m3 = h * derviY2(x + h / 2, y10[i] + m2 / 2, y20[i] + m2 / 2);
k4 = h * derviY1(x + h, y10[i] + k3, y20[i] + k3); m4 = h * derviY2(x + h, y10[i] + m3, y20[i] + m3);
y10[i + 1] = y10[i] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6; y20[i + 1] = y20[i] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6;
System.out.println("| " + x + " |" + " " + y10[i] + " " + "|" + " " + y20[i] + " " + "|"); } } }

Ответ

У вас h = 0. Почему? Ведь вы ясно написали, что h = (5-0)/10. А проблема вот в чём. Ваше выражение (5-0)/10 сначала преобразуется в int, а потом уже в double. То есть происходит всё как-то так (int)((5-0)/10). После этой операции мантисса отбрасывается и получается полноценный 0. Достаточно привести к double одну из переменных, как всё выражение будет обработано, как для переменных с плавающей точкой. К примеру так:
double h = ((double)b - a) / n;
В этих преобразованиях есть более интересные моменты, о которых хотелось бы рассказать, то это уже другая история.