Finite Element Approximation :: Galerkin Method

galerkin.cpp

/**
 * Galerkin Method
 * Solving for y''+ ay' + by + f(x) = 0 
 * a, b are constants
 * y(u1) = v1, y(u2) = v2 BVP 
 */

#include <iostream>
#include <cmath>
using namespace std;

struct node 
{
    int exp;
    float coeff;
    node *next;
    node (int e, float c) :
    exp(e), coeff(c), next(NULL) {}
};

node *insert(node *head, int exp, float coeff)
{
    if (!coeff) return head;
    if (head == NULL) return new node(exp, coeff);
    
    if (head->exp < exp) {
        node *newNode = new node(exp, coeff);
        newNode->next = head;
        return newNode;
    }

    node *iter = head, *prev;
    while (iter != NULL) {        
        if (iter->exp > exp && (iter->next == NULL || iter->next->exp < exp)) {
            node *newNode = new node(exp, coeff);
            newNode->next = iter->next;
            iter->next = newNode;
            return head;
        }
        if (iter->exp == exp) {
            iter->coeff += coeff;
            return head;
        }        
        prev = iter;
        iter = iter->next;
    }

    prev->next = new node(exp, coeff);

    return head;
}

node *add(node *poly1, node *poly2)
{
    node *newHead = NULL;

    if (poly1 == NULL) return poly2;
    if (poly2 == NULL) return poly1;

    node *iter1 = poly1, *iter2 = poly2, *iter;
    while (iter1 != NULL && iter2 != NULL) {
        if (iter1->exp > iter2-> exp) {
            if (newHead == NULL) {                
                iter = iter1;
                iter1 = iter1->next;    
                iter->next = NULL;
                newHead = iter;
            } else {
                iter->next = iter1;                
                iter1 = iter1->next;
                iter = iter->next;
                iter->next = NULL;
            }            
        } else if (iter1->exp < iter2-> exp) {
            if (newHead == NULL) {                
                iter = iter2;
                iter2 = iter2->next;
                iter->next = NULL;
                newHead = iter;
            } else {
                iter->next = iter2;
                iter2 = iter2->next;
                iter = iter->next;
                iter->next = NULL;
            }            
        } else {
            iter1->coeff += iter2->coeff;
            if (newHead == NULL) {       
                iter = iter1;
                iter1 = iter1->next;
                iter->next = NULL;
                newHead = iter;                
            } else {
                iter->next = iter1;
                iter1 = iter1->next;
                iter = iter->next;
                iter->next = NULL;
            }            
            node *temp = iter2;
            iter2 = iter2->next;            
            delete temp;
        }
    }

    if (iter1 != NULL) iter->next = iter1;
    if (iter2 != NULL) iter->next = iter2;

    return newHead;
}

node *multiply(node *poly1, node *poly2)
{
    if (poly1 == NULL || poly2 == NULL) return NULL;

    node *newHead = NULL;

    node *iter1 = poly1, *iter2, *iter;
    while (iter1 != NULL) {
        iter2 = poly2;
        while (iter2 != NULL) {
            newHead = insert(newHead, iter1->exp + iter2->exp,
                            (iter1->coeff) * (iter2->coeff));
            iter2 = iter2->next;
        }
        iter1 = iter1->next;
    }

    while (iter1 != NULL) {
        node *temp1 = iter1;
        iter1 = iter1->next;
        delete temp1;
    }

    while (iter2 != NULL) {
        node *temp2 = iter2;
        iter2 = iter2->next;
        delete temp2;
    }

    return newHead;
}

node *multiplyV(node *head, float x)
{   
    if (x == 0.0) return NULL;
    node *iter = head;
    while (iter != NULL) {
        iter->coeff *= x;
        iter = iter->next;
    }
    return head;
}

node *diff(node *head) 
{
    node *iter1 = head, *iter2 = NULL;
    while (iter1 != NULL) {
        if (iter1->exp == 0) {
            if (iter2 == NULL) return NULL;
            iter2->next = NULL;
            node *temp = iter1;
            delete temp;
            return head;
        } else {
            iter1->coeff *= iter1->exp;
            iter1->exp--;                        
        }
        iter2 = iter1;
        iter1 = iter1->next;
    }
    return head;
}

node *integ(node *head)
{
    node *iter1 = head;
    while (iter1 != NULL) {
        iter1->coeff /= (iter1->exp+1);
        iter1->exp++;
        iter1 = iter1->next;
    }
    return head;
}

float value(node *head, float x)
{
    float val = 0.0;
    while (head != NULL) {
        val += head->coeff * (pow(x, head->exp));
        head = head->next;
    }

    return val;
}

void printPoly(node *head)
{
    while (head != NULL) {
        cout << head->coeff;
        cout << "x" << head->exp << ((head->next != NULL) ? " + " : ""); 
        head = head->next;
    }
    cout << endl;
}

node *copy(node *head)
{
    node *newHead = NULL;
    while (head != NULL) {
        newHead = insert(newHead, head->exp, head->coeff);
        head = head->next;
    }
    return newHead;
}

int main (int argc, char *argv[])
{
    float a = 0, b = 1;
    
    node *f = NULL;
    f = insert(f, 1, 2);
    f = insert(f, 2, -2);

    float u1, u2, v1, v2;
    u1 = 0, v1 = 0;
    u2 = 1, v2 = 0.5;

    node *phi[2] = {NULL};
    phi[0] = insert(phi[0], 0, v1 - u1 * ((v2-v1)*1.0/(u2-u1)));
    phi[0] = insert(phi[0], 1, (v2-v1)*1.0/(u2-u1));

    node *p1 = NULL, *p2 = NULL;
    p1 = insert(p1, 1, 1);
    p2 = insert(p2, 0, -u1);
    p2 = insert(p2, 1, 1);
    p2 = insert(p2, 0, -u2);
    phi[1] = multiply(p1, p2);

    node *dphi[2] = {NULL};
    dphi[0] = diff(copy(phi[0]));
    dphi[1] = diff(copy(phi[1]));

    node *ddphi[2] = {NULL};
    ddphi[0] = diff(copy(dphi[0]));
    ddphi[1] = diff(copy(dphi[1]));

    node *iphi[2] = {NULL};    
    iphi[0] = 
    add(
        add(
            integ(multiply(copy(ddphi[0]), copy(phi[1]))),
            multiplyV(integ(multiply(copy(dphi[0]), copy(phi[1]))), a)
        ),
        add(
            multiplyV(integ(multiply(copy(phi[0]), copy(phi[1]))), b),
            integ(multiply(copy(f), copy(phi[1])))
        )
    );

    iphi[1] = 
    add(
        add(
            integ(multiply(copy(ddphi[1]), copy(phi[1]))),
            multiplyV((integ(multiply(copy(dphi[1]), copy(phi[1])))), a)
        ),
        multiplyV((integ(multiply(copy(phi[1]), copy(phi[1])))), b)
    );

    float lhs = value(iphi[0], u2)-value(iphi[0], u1),
          rhs = value(iphi[1], u2)-value(iphi[1], u1);

    float C = rhs/lhs;
    
    node *ans = add(copy(phi[0]), multiplyV(copy(phi[1]), C));
    printPoly(ans);
    
    return 0;
}