Как найти корни кубического уравнения в питоне

Here, I put a cubic equation (with complex coefficients) solver.

#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>

using namespace std;

#define PI 3.141592

long double complex_multiply_r(long double xr, long double xi, long double yr, long double yi) {
    return (xr * yr - xi * yi);
}

long double complex_multiply_i(long double xr, long double xi, long double yr, long double yi) {
    return (xr * yi + xi * yr);
}

long double complex_triple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
    return (xr * yr * zr - xi * yi * zr - xr * yi * zi - xi * yr * zi);
}

long double complex_triple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi) {
    return (xr * yr * zi - xi * yi * zi + xr * yi * zr + xi * yr * zr);
}

long double complex_quadraple_multiply_r(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
    long double z1r, z1i, z2r, z2i;    
    z1r = complex_multiply_r(xr, xi, yr, yi);
    z1i = complex_multiply_i(xr, xi, yr, yi);
    z2r = complex_multiply_r(zr, zi, wr, wi);
    z2i = complex_multiply_i(zr, zi, wr, wi);
    return (complex_multiply_r(z1r, z1i, z2r, z2i));
}

long double complex_quadraple_multiply_i(long double xr, long double xi, long double yr, long double yi, long double zr, long double zi, long double wr, long double wi) {
    long double z1r, z1i, z2r, z2i;
    z1r = complex_multiply_r(xr, xi, yr, yi);
    z1i = complex_multiply_i(xr, xi, yr, yi);
    z2r = complex_multiply_r(zr, zi, wr, wi);
    z2i = complex_multiply_i(zr, zi, wr, wi);
    return (complex_multiply_i(z1r, z1i, z2r, z2i));
}

long double complex_divide_r(long double xr, long double xi, long double yr, long double yi) {
    return ((xr * yr + xi * yi) / (yr * yr + yi * yi));
}

long double complex_divide_i(long double xr, long double xi, long double yr, long double yi) {
    return ((-xr * yi + xi * yr) / (yr * yr + yi * yi));
}

long double complex_root_r(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (sqrt(r) * cos(theta / 2.0));
    }
    else {
        return 0.0;
    }

}    

long double complex_root_i(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (sqrt(r) * sin(theta / 2.0));
    }
    else {
        return 0.0;
    }
}    

long double complex_cuberoot_r(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (pow(r, 1.0 / 3.0) * cos(theta / 3.0));
    }
    else {
        return 0.0;
    }
}    

long double complex_cuberoot_i(long double xr, long double xi) {
    long double r, theta;
    r = sqrt(xr*xr + xi*xi);
    if (r != 0.0) {
        if (xr >= 0 && xi >= 0) {
            theta = atan(xi / xr);
        }
        else if (xr < 0 && xi >= 0) {
            theta = PI - abs(atan(xi / xr));
        }
        else if (xr < 0 && xi < 0) {
            theta = PI + abs(atan(xi / xr));
        }
        else {
            theta = 2.0 * PI + atan(xi / xr);
        }
        return (pow(r, 1.0 / 3.0) * sin(theta / 3.0));
    }
    else {
        return 0.0;
    }
}    

void main() {
    long double a[2], b[2], c[2], d[2], minusd[2];
    long double r, theta;
    cout << "ar?";
    cin >> a[0];
    cout << "ai?";
    cin >> a[1];
    cout << "br?";
    cin >> b[0];
    cout << "bi?";
    cin >> b[1];
    cout << "cr?";
    cin >> c[0];
    cout << "ci?";
    cin >> c[1];
    cout << "dr?";
    cin >> d[0];
    cout << "di?";
    cin >> d[1];

    if (b[0] == 0.0 && b[1] == 0.0 && c[0] == 0.0 && c[1] == 0.0) {
        if (d[0] == 0.0 && d[1] == 0.0) {
            cout << "x1r: 0.0 n";
            cout << "x1i: 0.0 n";
            cout << "x2r: 0.0 n";
            cout << "x2i: 0.0 n";
            cout << "x3r: 0.0 n";
            cout << "x3i: 0.0 n";
        }
        else {
                minusd[0] = -d[0];
                minusd[1] = -d[1];
                r = sqrt(minusd[0]*minusd[0] + minusd[1]*minusd[1]);
                if (minusd[0] >= 0 && minusd[1] >= 0) {
                    theta = atan(minusd[1] / minusd[0]);
                }
                else if (minusd[0] < 0 && minusd[1] >= 0) {
                    theta = PI - abs(atan(minusd[1] / minusd[0]));
                }
                else if (minusd[0] < 0 && minusd[1] < 0) {
                    theta = PI + abs(atan(minusd[1] / minusd[0]));
                }
                else {
                    theta = 2.0 * PI + atan(minusd[1] / minusd[0]);
                }
                cout << "x1r: " << pow(r, 1.0 / 3.0) * cos(theta / 3.0) << "n";
                cout << "x1i: " << pow(r, 1.0 / 3.0) * sin(theta / 3.0) << "n";
                cout << "x2r: " << pow(r, 1.0 / 3.0) * cos((theta + 2.0 * PI) / 3.0) << "n";
                cout << "x2i: " << pow(r, 1.0 / 3.0) * sin((theta + 2.0 * PI) / 3.0) << "n";
                cout << "x3r: " << pow(r, 1.0 / 3.0) * cos((theta + 4.0 * PI) / 3.0) << "n";
                cout << "x3i: " << pow(r, 1.0 / 3.0) * sin((theta + 4.0 * PI) / 3.0) << "n";
            }
        }
        else {
        // find eigenvalues
        long double term0[2], term1[2], term2[2], term3[2], term3buf[2];
        long double first[2], second[2], second2[2], third[2];
        term0[0] = -4.0 * complex_quadraple_multiply_r(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
        term0[1] = -4.0 * complex_quadraple_multiply_i(a[0], a[1], c[0], c[1], c[0], c[1], c[0], c[1]);
        term0[0] += complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
        term0[1] += complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], c[0], c[1], c[0], c[1]);
        term0[0] += -4.0 * complex_quadraple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
        term0[1] += -4.0 * complex_quadraple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1], d[0], d[1]);
        term0[0] += 18.0 * complex_quadraple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
        term0[1] += 18.0 * complex_quadraple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1], d[0], d[1]);
        term0[0] += -27.0 * complex_quadraple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
        term0[1] += -27.0 * complex_quadraple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1], d[0], d[1]);
        term1[0] = -27.0 * complex_triple_multiply_r(a[0], a[1], a[0], a[1], d[0], d[1]);
        term1[1] = -27.0 * complex_triple_multiply_i(a[0], a[1], a[0], a[1], d[0], d[1]);
        term1[0] += 9.0 * complex_triple_multiply_r(a[0], a[1], b[0], b[1], c[0], c[1]);
        term1[1] += 9.0 * complex_triple_multiply_i(a[0], a[1], b[0], b[1], c[0], c[1]);
        term1[0] -= 2.0 * complex_triple_multiply_r(b[0], b[1], b[0], b[1], b[0], b[1]);
        term1[1] -= 2.0 * complex_triple_multiply_i(b[0], b[1], b[0], b[1], b[0], b[1]);
        term2[0] = 3.0 * complex_multiply_r(a[0], a[1], c[0], c[1]);
        term2[1] = 3.0 * complex_multiply_i(a[0], a[1], c[0], c[1]);
        term2[0] -= complex_multiply_r(b[0], b[1], b[0], b[1]);
        term2[1] -= complex_multiply_i(b[0], b[1], b[0], b[1]);
        term3[0] = complex_multiply_r(term1[0], term1[1], term1[0], term1[1]);
        term3[1] = complex_multiply_i(term1[0], term1[1], term1[0], term1[1]);
        term3[0] += 4.0 * complex_triple_multiply_r(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
        term3[1] += 4.0 * complex_triple_multiply_i(term2[0], term2[1], term2[0], term2[1], term2[0], term2[1]);
        term3buf[0] = term3[0];
        term3buf[1] = term3[1];
        term3[0] = complex_root_r(term3buf[0], term3buf[1]);
        term3[1] = complex_root_i(term3buf[0], term3buf[1]);

        if (term0[0] == 0.0 && term0[1] == 0.0 && term1[0] == 0.0 && term1[1] == 0.0) {
            cout << "x1r: " << -pow(d[0], 1.0 / 3.0) << "n";
            cout << "x1i: " << 0.0 << "n";
            cout << "x2r: " << -pow(d[0], 1.0 / 3.0) << "n";
            cout << "x2i: " << 0.0 << "n";
            cout << "x3r: " << -pow(d[0], 1.0 / 3.0) << "n";
            cout << "x3i: " << 0.0 << "n";
        }
        else {
            // eigenvalue1
            first[0] = complex_divide_r(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1]), 3.0 * pow(2.0, 1.0 / 3.0) * a[0], 3.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(pow(2.0, 1.0 / 3.0) * term2[0], pow(2.0, 1.0 / 3.0) * term2[1], 3.0 * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x1r: " << first[0] - second[0] - third[0] << "n";
            cout << "x1i: " << first[1] - second[1] - third[1] << "n";

            // eigenvalue2
            first[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, -sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x2r: " << -first[0] + second[0] - third[0] << "n";
            cout << "x2i: " << -first[1] + second[1] - third[1] << "n";

            // eigenvalue3
            first[0] = complex_divide_r(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            first[1] = complex_divide_i(complex_multiply_r(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), complex_multiply_i(1.0, sqrt(3.0), complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 6.0 * pow(2.0, 1.0 / 3.0) * a[0], 6.0 * pow(2.0, 1.0 / 3.0) * a[1]);
            second[0] = complex_divide_r(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            second[1] = complex_divide_i(complex_multiply_r(1.0, -sqrt(3.0), term2[0], term2[1]), complex_multiply_i(1.0, -sqrt(3.0), term2[0], term2[1]), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_r(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])), 3.0 * pow(2.0, 2.0 / 3.0) * complex_multiply_i(a[0], a[1], complex_cuberoot_r(term3[0] + term1[0], term3[1] + term1[1]), complex_cuberoot_i(term3[0] + term1[0], term3[1] + term1[1])));
            third[0] = complex_divide_r(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            third[1] = complex_divide_i(b[0], b[1], 3.0 * a[0], 3.0 * a[1]);
            cout << "x3r: " << -first[0] + second[0] - third[0] << "n";
            cout << "x3i: " << -first[1] + second[1] - third[1] << "n";
        }
    }

    int end;
    cin >> end;
}

Источник
for (int x = -100; x

You can’t perform that action at this time.

You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.

Библиотека Sympy: символьные вычисления в Python

Что такое SymPy ? Это библиотека символьной математики языка Python. Она является реальной альтернативой таким математическим пакетам как Mathematica или Maple и обладает очень простым и легко расширяемым кодом. SymPy написана исключительно на языке Python и не требует никаких сторонних библиотек.

Документацию и исходный код этой библиотеки можно найти на ее официальной странице.

Первые шаги с SymPy

Используем SymPy как обычный калькулятор

В библиотеке SymPy есть три встроенных численных типа данных: Real , Rational и Integer . С Real и Integer все понятно, а класс Rational представляет рациональное число как пару чисел: числитель и знаменатель рациональной дроби. Таким образом, Rational(1, 2) представляет собой 1/2 , а, например, Rational(5, 2) — соответственно 5/2 .

Библиотека SymPy использует библиотеку mpmath , что позволяет производить вычисления с произвольной точностью. Таким образом, ряд констант (например, пи, e), которые в данной библиотеке рассматриваются как символы, могут быть вычислены с любой точностью.

Как можно заметить, функция evalf() дает на выходе число с плавающей точкой.

В SymPy есть также класс, представляющий такое понятие в математике, как бесконечность. Он обозначается следующим образом: oo .

Символы

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

После их задания, с ними можно производить различные манипуляции.

С символами можно производить преобразования с использованием некоторых операторов языка Python. А именно, арифметических ( + , -` , «* , ** ) и логических ( & , | ,

Библиотека SymPy позволяет задавать форму вывода результатов на экран. Обычно мы используем формат такого вида:

Алгебраические преобразования

SymPy способна на сложные алгебраические преобразования. Здесь мы рассмотрим наиболее востребованные из них, а именно раскрытие скобок и упрощение выражений.

Раскрытие скобок

Чтобы раскрыть скобки в алгебраических выражениях, используйте следующий синтаксис:

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

Упрощение выражений

Если вы хотите привести выражение к более простому виду (возможно, сократить какие-то члены), то используйте функцию simplify .

Также надо сказать, что для определенных видов математических функций существуют альтернативные, более конкретные функции для упрощения выражений. Так, для упрощения степенных функций есть функция powsimp , для тригонометрических — trigsimp , а для логарифмических — logcombine , radsimp .

Вычисления

Вычисления пределов

Для вычисления пределов в SymPy предусмотрен очень простой синтаксис, а именно limit(function, variable, point) . Например, если вы хотите вычислить предел функции f(x) , где x -> 0 , то надо написать limit(f(x), x, 0) .

Также можно вычислять пределы, которые стремятся к бесконечности.

Дифференцирование

Для дифференцирования выражений в SymPy есть функция diff(func, var) . Ниже даны примеры ее работы.

Проверим результат последней функции при помощи определения производной через предел.

tan 2 (𝑥)+1 Результат тот же.

Также при помощи этой же функции могут быть вычислены производные более высоких порядков. Синтаксис функции будет следующим: diff(func, var, n) . Ниже приведено несколько примеров.

Разложение в ряд

Для разложения выражения в ряд Тейлора используется следующий синтаксис: series(expr, var) .

Интегрирование

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

Также несложно посчитать интеграл и от специальных функций. Возьмем, например, функцию Гаусса:

Результат вычисления можете посмотреть сами. Вот примеры вычисления определенных интегралов.

Также можно вычислять определенные интегралы с бесконечными пределами интегрирования (несобственные интегралы).

Решение уравнений

При помощи SymPy можно решать алгебраические уравнения с одной или несколькими переменными. Для этого используется функция solveset() .

Как можно заметить, первое выражение функции solveset() приравнивается к 0 и решается относительно х . Также возможно решать некоторые уравнения с трансцендентными функциями.

Системы линейных уравнений

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

Факторизация

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

Булевы уравнения

Также в SymPy реализована возможность решения булевых уравнений, что по сути означает проверку булевого выражения на истинность. Для этого используется функция satisfiable() .

Данный результат говорит нам о том, что выражение (x & y) будет истинным тогда и только тогда, когда x и y истинны. Если выражение не может быть истинным ни при каких значениях переменных, то функция вернет результат False .

Линейная алгебра

Матрицы

Матрицы в SymPy создаются как экземпляры класса Matrix :

В отличие от NumPy , мы можем использовать в матрицах символьные переменные:

И производить с ними разные манипуляции:

Дифференциальные уравнения

При помощи библиотеки SymPy можно решать некоторые обыкновенные дифференциальные уравнения. Для этого используется функция dsolve() . Для начала нам надо задать неопределенную функцию. Это можно сделать, передав параметр cls=Function в функцию symbols() .

Теперь f и g заданы как неопределенные функции. мы можем в этом убедиться, просто вызвав f(x) .

Теперь решим следующее дифференциальное уравнение:

Чтобы улучшить решаемость и помочь этой функции в поиске решения, можно передавать в нее определенные ключевые аргументы. Например, если мы видим, что это уравнение с разделяемыми переменными, то мы можем передать в функцию аргумент hint=’separable’ .

Бесплатные кодинг марафоны с ревью кода

Наш телеграм канал проводит бесплатные марафоны по написанию кода на Python с ревью кода от преподавателя

источники:

http://gist.github.com/LygutaKsusha/0f76f728bdbbf074be99

Библиотека Sympy: символьные вычисления в Python

Источник

You could implement the cubic formula

enter image description here

this Youtube video from mathologer could help understand it.

Based on that, the cubic function for ax^3 + bx^2 + cx + d = 0 can be written like this:

def cubic(a,b,c,d):
    n = -b**3/27/a**3 + b*c/6/a**2 - d/2/a
    s = (n**2 + (c/3/a - b**2/9/a**2)**3)**0.5
    r0 = (n-s)**(1/3)+(n+s)**(1/3) - b/3/a
    r1 = (n+s)**(1/3)+(n+s)**(1/3) - b/3/a
    r2 = (n-s)**(1/3)+(n-s)**(1/3) - b/3/a
    return (r0,r1,r2)

The simplified version of the formula only needs to get c and d as parameters (aka p and q) and can be implemented like this:

def cubic(p,q):
    n = -q/2
    s = (q*q/4+p**3/27)**0.5
    r0 = (n-s)**(1/3)+(n+s)**(1/3)
    r1 = (n+s)**(1/3)+(n+s)**(1/3)
    r2 = (n-s)**(1/3)+(n-s)**(1/3)
    return (r0,r1,r2)

print(cubic(-15,-126))
(5.999999999999999, 9.999999999999998, 2.0)

I’ll let you mix in complex number operations to properly get all 3 roots

Источник

0 / 0 / 0

Регистрация: 21.09.2021

Сообщений: 1

1

Решение кубического уравнения

21.09.2021, 20:59. Показов 10932. Ответов 9

Студворк — интернет-сервис помощи студентам

Ребята можете, пожалуйста, написать рабочий простой код, который решает кубическое уравнение. На ввод подаются коэффициенты a, b, c ,d. При вводе этих коэффициентов код должен решить и выдать корни решения уравнения (комплексные числа допускаются). Кратко: нужен код, который выдаст корни кубического уравнения при любых коэффициентах ввода.



0

enx

1182 / 758 / 277

Регистрация: 05.09.2021

Сообщений: 1,772

21.09.2021, 21:05

2

Проверяйте,

Python 
1
2
3
4
5

from sympy import solve, Symbol
 
x = Symbol('x')
a, b, c, d = (int(input()) for i in 'abcd')
print(solve(a * x ** 3 + b * x ** 2 + c * x + d, x))



1

Gdez

Эксперт Python

7258 / 4047 / 1780

Регистрация: 27.03.2020

Сообщений: 6,875

21.09.2021, 21:18

3

_Nayt_,

Python 
1
2
3
4

from numpy.polynomial import Polynomial
a,b,c,d = map(float(input()
x = Polynomial([d,c,b,a])
print(x.roots())

a -> коэффициент при x^3

Добавлено через 1 минуту
Почти любой полином



1

1182 / 758 / 277

Регистрация: 05.09.2021

Сообщений: 1,772

21.09.2021, 21:21

4

Цитата
Сообщение от Gdez
Посмотреть сообщение

float(

Вот за это не люблю я map



0

Эксперт Python

7258 / 4047 / 1780

Регистрация: 27.03.2020

Сообщений: 6,875

21.09.2021, 21:27

5

enx, — ночерь

Добавлено через 3 минуты
enx, вообще: насколько мап медленнее генератора?



0

enx

21.09.2021, 21:35

Не по теме:

Цитата
Сообщение от Gdez
Посмотреть сообщение

медленнее генератора

Я не проверял никогда, народ говорит что примерно на 1/3.



0

Эксперт Python

4607 / 2028 / 359

Регистрация: 17.03.2012

Сообщений: 10,086

Записей в блоге: 6

22.09.2021, 08:57

7

Цитата
Сообщение от Gdez
Посмотреть сообщение

насколько мап медленнее генератора?

Так map же и есть генератор (в 3-м питоне)? Или вы обо что?



0

Catstail

Модератор

Эксперт функциональных языков программированияЭксперт Python

35427 / 19452 / 4071

Регистрация: 12.02.2012

Сообщений: 32,488

Записей в блоге: 13

22.09.2021, 10:35

8

Лучший ответ Сообщение было отмечено Fudthhh как решение

Решение

enx, Gdez, вам, увы — незачет… А вот Ондржею Чертику и Аарону Мереру (это авторы sympy) — разумеется 5++

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

Python 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

def sign(x):  # вроде бы есть в numpy, но и написать несложно...
    if x>0:
        return 1
    elif x<0:
        return -1
    else:
        return 0
        
def dicho(f,a,b,eps=1.0e-14):  # Дихотомия. Один вещественный корень у кубического полинома 
    fa=f(a)                               # есть всегда. Его и найдем...
    fb=f(b)
    while True:
        c=0.5*(a+b)
        if abs(b-a)<eps:
            return c
        fc=f(c)
        if abs(fc)<=eps:
            return c
        if sign(fa)*sign(fc)<0:
            b=c
            fb=fc
        else:
            a=c
            fa=fc
            
def div_poly(p,a):  # делим исходный полином на (x-a) и получаем квадратный
    r=[0,0,0]          # трехчлен, поставляющий два недостающих корня
    r[2]=p[3]
    r[1]=p[2]+a*p[3]
    r[0]=(p[1]+a*(p[2]+a*p[3]))
    return tuple(r)
    
def solve_qube(p):  # парадная функция
    q=max(p)
    left=-abs(q)/abs(p[3])  # границы вещественного корня
    right=-left
    x1=dicho(lambda x: p[3]*x**3+p[2]*x**2+p[1]*x+p[0],left,right)
    (c,b,a)=div_poly(p,x1)
    d=b**2-4*a*c
    x2=(-b+d**0.5)/(2*a)
    x3=(-b-d**0.5)/(2*a)
    return (x1,x2,x3)
    
print(solve_qube([12,4,8,5]))

Вывод:

(-1.8623879412594293, (0.13119397062971466+1.1275887015236796j), (0.13119397062971452-1.1275887015236796j))



2

1182 / 758 / 277

Регистрация: 05.09.2021

Сообщений: 1,772

22.09.2021, 10:49

9

Цитата
Сообщение от Catstail
Посмотреть сообщение

вам, увы — незачет

Да вы сразу байткод пишите, что уж тут )))



1

Модератор

Эксперт функциональных языков программированияЭксперт Python

35427 / 19452 / 4071

Регистрация: 12.02.2012

Сообщений: 32,488

Записей в блоге: 13

22.09.2021, 10:54

10

enx, а смысл? Это рутинная процедура, ничего не дающая уму. Как, впрочем, и вызов стандартных функций, написанных другими. А мой подход полезнее для обучаемого.



0
Источник


This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters

Show hidden characters

import java.util.*;
import java.io.*;
/**
* Created by Ксю on 07.09.2014.
*/
public class Kub {
public static void main(String[] args) throws IOException {
// определяем значения коефициентов
int a1=1; int b1 =-3; int c1= 0; int d1 =0;
int a2=3; int b2=-15; int c2=18; int d2=0;
int a3 =1; int b3 =-7; int c3 =-33; int d3 =135;
// передаем в метод значения
System.out.println(«Kub1 is: «);
kub(a1,b1,c1,d1);
System.out.println(«Kub2 is: «);
kub(a2,b2,c2,d2);
System.out.println(«Kub3 is: «);
kub(a3,b3,c3,d3);
}
// метод для решения кубического уравнения
public static void kub(int a, int b, int c, int d) {
for (int x = -100; x <= 100; x++) {
double result = a * Math.pow(x, 3) + b * x * x + c * x + d;
if (result == 0){
System.out.println(«X = » + x);
}
}
}
}
Источник

Понравилась статья? Поделить с друзьями:
  • Как составить план работы службы
  • Как найти показания счетчика электроэнергии
  • Fltmgr sys синий экран windows 10 x64 как исправить
  • Как составить штатное расписание с почасовой оплатой труда
  • Как найти разность четырех чисел