lab3.cpp

The decision of a boundary problem a method of Galerkin

// lab3.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <math.h>
#include <stdio.h>

const double a = 0.0;
const double b = 2.;
const double eps = 0.0001;

double p(double x)
{
    return exp(3 - ( 3 + x * x + 2 * x ) * exp(-x));
}

double q(double x)
{
    return - (-exp(-x)/( 1 + x*x ) * p(x));
}

double f(double x)
{
    return - (15 * (x*x + 1) * p(x));
}


double Phi(double x)
{
    if (fabs(x) <= 1) return (1 - fabs(x));
    else return 0;
}

double PhiK(double x, double xk, double h)
{
    return Phi((x - xk)/h);
}

double DPhi(double x)
{
    if (-1 <= x && x <=0) return 1;
    if (0 <= x && x <= 1) return -1;
    return 0;
}

double DPhiK(double x, double xk, double h)
{
    return DPhi((x - xk)/h)/h;
}

double F(double x, int k1, int k2, double h)
{
    return p(x)*DPhiK(x,a+k1*h,h)*DPhiK(x,a+k2*h,h) + q(x)*PhiK(x,a+k1*h,h)*PhiK(x,a+k2*h,h);
}

double Gauss (int k1, int k2, double h)
{
    double c[4], t[4], sum=0, hm;
    int i, k, ng, mg;
    
    mg=1000;
    ng=4;
    c[0] = c[3] = 0.34785484513774368;
    c[1] = c[2] = 0.6521451548625432;
    
    t[0] = -(t[3] = 0.8611363115940492);
    t[1] = -(t[2] = 0.3399810435848646);
    
    hm = (b-a)/mg;
    
    for ( k = 1; k <= mg; k++ )
        for ( i = 1; i <= ng; i++)
            sum+= c[i-1] * F(a+(2*k-1)*hm/2+hm/2*t[i-1],k1,k2,h);
        
        sum*= hm/2;
        return sum;
}

double FPhi(double x, int k, double h)
{
    return PhiK(x,a+k*h,h)*f(x);
}

double Scal(int l, double h)
{
    double c[4], t[4], sum=0, hm;
    int i, k, ng, mg;
    
    mg=1000;
    ng=4;
    c[0] = c[3] = 0.34785484513774368;
    c[1] = c[2] = 0.6521451548625432;
    
    t[0] = -(t[3] = 0.8611363115940492);
    t[1] = -(t[2] = 0.3399810435848646);
    
    hm = (b-a)/mg;
    
    for ( k = 1; k <= mg; k++ )
        for ( i = 1; i <= ng; i++)
            sum+= c[i-1] * FPhi( a + ( 2 * k - 1 ) * hm / 2 + hm / 2 * t [ i - 1 ], l, h );
        
    sum *= hm/2;
    return sum;
}



void main()
{
    double  h,d, k0, k1, n0, n1, 
        m0 = 0.0, // m0
        m1 = 0.0, // m1
        M0 = 1.0, // 斐0
        M1 = 0.0; // 吵1
    int i, n, N = 20, k, z=0, j;
    
    double *Y1 = new double[N+1];
    double *Y2 = new double[N+1];
    double *NV = new double[N+1];
    
    for (i = 0; i <= N; i++)
    {
        Y1[i] = 10;
        Y2[i] = 1000;
    }
    
    k = 5;
    
    do
    {
        z=0;
        n = k*N;
        
        double *A = new double[n+1];
        double *B = new double[n+1];
        double *C = new double[n+1];
        double *D = new double[n+1];
        
        double *ksi = new double[n+1];
        
        double *Alpha = new double[n+1];
        double *Betha = new double[n+1];
        
        h = (b-a)/n;
        
        
        B[0] = Gauss(0,0,h);
        C[0] = Gauss(1,0,h);
        for(i = 1; i <= n; i++)
        {
            A[i] = Gauss(i-1,i,h);
            B[i] = Gauss(i,i,h);
            C[i] = Gauss(i+1,i,h);
            D[i] = Scal(i, h);
        }
        D[0] = m0*p(a) + Scal(0,h);
        D[n] = m1*p(b) + Scal(n,h);
        
        k0 = -C[0]/(B[0] + M0 * p(a));
        k1 = -A[n]/(B[n] + M1*p(b));
        n0 = D[0]/(B[0] + M0*p(a));
        n1 = D[n]/(B[n] + M1*p(b));
        
        Alpha[0] = k0;
        Betha[0] = n0;
        
        for(i = 1; i <= n-1; i++)
        {
            Alpha[i] = (-C[i]) / (B[i] + A[i]*Alpha[i-1]);
            Betha[i] = (D[i] - A[i]*Betha[i-1]) / (B[i] + A[i]*Alpha[i-1]);
        }
        
        ksi[n] = m1;
        for(i = n-1; i >= 0; i--)
        {
            ksi[i] = Alpha[i]*ksi[i+1] + Betha[i];
        }
        
        for (i = 0; i <= N; i++)
        {
            d=0;
            for (j=0;j<=n;j++)
                d+=ksi[j]*PhiK(a+i*k*h,a+j*h,h);
            Y1[i] = d;
        }
        
        for (i = 0; i <= N; i++)
        {
            NV[ i ] = fabs(Y1[i] - Y2[i]);
            if ( NV[i] < eps )
                z++;
        }
            
        for (i = 0; i <= N; i++)
        {
            Y2[i] = Y1[i];
        }
        
        delete [] A;
        delete [] B;
        delete [] C;
        delete [] D;
        
        delete [] ksi;
        delete [] Alpha;
        delete [] Betha;
        
        k*=2;
    } while (z < 3);
    
 
    k/=2;
    printf(" i |  X[i]  |     Y[i]     |   Nevjazka |\n");
    printf("-----------------------------------------\n");
    for (i=0;i<=N;i++)
        printf("%-2d | %-6.4f | %-+11.9f | %-10.8f |\n", i, a + k * i * h, Y2[i], NV[i]);//cout << a+k*i*h << "   "<<Y2[i]<<endl;
    
    delete [] Y1;
    delete [] Y2;
    delete [] NV;

}