Разложение синуса/косинуса в ряд по Тейлору

--[[Участник:Artix|Artix]] 11:40, 26 января 2011 (UTC)

Классическая задача матанализа - разложение функции в ряд. 

Кроме sin(x) и cos(x) я попытался разложить еще и другие функции - 
и тут наткнулся на сложности в реализации. Для вычисления тангенса понадобилось
реализовать биноминальные коэффициенты Ньютона и числа Бернулли, а для вычисления
секанса - еще и числа Эйлера. С последними долго не получалось, т.к. пришлось всё
думать самому по формулам. В интернете никаких вычислительных примеров на эту тему
я не нашел. Суть сложностей состоит в том - что я не использовал модуль math.h для
вычислений, а пользовался исключительно обычными числами и дробями, без применения
логарифмов и т.п.

/*///////////////////////////////////////////////////////////////

 Разложение sin(), cos() по Тейлору, GNU C++, 2011-01-01
 
 Artix, master{@}7masterov.ru, skype:artixmaster
 
 * Error in code? Nothing is perfect!
 * Free source for free Linux, use it for free!
 * Please, do not remove this comment!

///////////////////////////////////////////////////////////////*/

#include <iostream>
#include <math.h>

using namespace std;

const long N=12;


// Модуль
double mabs(double x)
{
    return (x>0)?x:-x;
}

// Степерь x^y
double mpow(double x, long y)
{
    double r = 1;
    while((y--)>0) r*=x;
    return r;
}

// Факториал x!
long fact(long x)
{
    if (x<=1) return 1;
    else return x*fact(x-1);
}

// Биноминальный коэффициент Ньютона C(n/k)
double bink(long n, long k)
{
    return 1.*fact(n)/fact(k)/fact(n-k);
}

// Числа Бернулли B(n)
double bernoulli(long n)
{
    if (n<=0) return 1;
    else {
        double s = 0;
        for(long k=1;k<=n;k++) s+=bink(n+1,k+1)*bernoulli(n-k);
        return -1./(n+1)*s;
    }
}

// Ряд синуса
double seq_sin(double x)
{
    double r = 0;
    for(long n=0;n<N;n++) {
        r+=mpow(-1,n)*mpow(x,2*n+1)/fact(2*n+1);
    }
    return r;
}

// Ряд косинуса
double seq_cos(double x)
{
    double r = 0;
    for(long n=0;n<N;n++) {
        r+=mpow(-1,n)*mpow(x,2*n)/fact(2*n);
    }
    return r;
}

// Ряд тангенса
double seq_tan(double x)
{
    double r = 0;
    for(long n=1;n<N;n++) {
        r+=mpow(2,2*n)*(mpow(2,2*n)-1)*mabs(bernoulli(2*n))/fact(2*n)*mpow(x,2*n-1);
    }
    return r;
}

// Ряд котангенса
double seq_ctg(double x)
{
    double r = 1./x;
    for(long n=1;n<=N;n++) {
        r+=mpow(-1,n)*mpow(2,2*n)*(bernoulli(2*n))*mpow(x,2*n-1)/fact(2*n);
    }
    return r;
}

// Числа Эйлера 1 рода - не используются
double euler1(long n, long k)
{
    double res=0;
    for(long j=0;j<=k;j++)
        res+=bink(n+1,j)*mpow(1+k-j,n)*mpow(-1,j);
    return res;
}

// Числа Эйлера
double euler(long n,int flag=1)
{
    if (n==0) return 1;
    if (flag) {
        if (n%2==1) return 0;
        n/=2;
    }
    double res=0;
    for(long j=0;j<=n-1;j++)
        res+=bink(2*n,2*j)*euler(j,0);
    return -res;
}

// Ряд секанса
double seq_sec(double x)
{
    double r = 1;
    for(long n=1;n<N;n++) {
        r+=mabs(euler(2*n))/fact(2*n)*mpow(x,2*n);
    }
    return r;
}

int main()
{
    cout << "sin(pi/5)=" << sin(M_PI/5) << endl;
    cout << "seq_sin(pi/5)="<< seq_sin(M_PI/5) << endl;

    cout << "cos(pi/5)=" << cos(M_PI/5) << endl;
    cout << "seq_cos(pi/5)="<< seq_cos(M_PI/5) << endl;

    cout << "tan(pi/5)=" << tan(M_PI/5) << endl;
    cout << "seq_tan(pi/5)="<< seq_tan(M_PI/5) << endl;

    cout << "ctg(pi/5)=" << 1./tan(M_PI/5) << endl;
    cout << "seq_ctg(pi/5)="<< seq_ctg(M_PI/5) << endl;

    cout << "sec(pi/5)=" << 1./cos(M_PI/5) << endl;
    cout << "seq_sec(pi/5)="<< seq_sec(M_PI/5) << endl;
    return 1;
}
// eof