Laboratory Works

Laboratory work No. 5

Quadratic

(Python code of the problem)

The problem was solved using LAGRANGE priciple. The values of x0, x1 and x2 as the values of y0, y1 and y2 are changed for points a), b) and c). The images are shown below the code, for each of a), b) and c).

import matplotlib.pyplot as plt
import numpy as np

pi = 3.14

_x0 = 3
_y0 = 1.0986123

_x1 = 4
_y1 = 1.3862944

_x2 = 5
_y2 = 1.6094379

_x = np.arange(-500, 500, 0.1);



def L_function(x, x1, x2, x0, y0, y1, y2):
    L20 = ((x - x1) * (x - x2)) / ((x0 - x1) * (x0 - x2))
    L21 = ((x - x0) * (x - x2)) / ((x1 - x0) * (x1 - x2))
    L22 = ((x - x0) * (x - x1)) / ((x2 - x0) * (x2 - x1))

    L = y0 * L20 + y1 * L21 + y2 * L22
    return L

_y = L_function(_x, _x1, _x2, _x0, _y0, _y1, _y2)
print (_y)


plt.plot(_y)
plt.show()

Results:

a)

b)

c)

Using Lagrangian interpolation, we were allowe to simulate the graph of the function. Even if we obtained some results, Lagrange interpolation is essentially NEVER a good choice for interpolation, because in reality, the obrained graph cand be far from real function graph. Therefore, is a good way to INTRODUCE ideas of interpolation, and sometimes to prove some simple results. But for precision you need something more.

Cubic

(Python code of the problem)

Same principle as in exercise Nr. 1

import matplotlib.pyplot as plt
import numpy as np

pi = 3.14

_x0 = 1
_y0 = 0

_x1 = 3
_y1 = 1.0986123

_x2 = 4
_y2 = 1.3862944

_x3 = 6
_y3 = 1.7917595

_x = np.arange(-50, 50, 0.1);



def L_function(x, x1, x2, x0, x3, y0, y1, y2, y3):
    L30 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3))
    L31 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3))
    L32 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3))
    L33 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2))

    L = y0 * L30 + y1 * L31 + y2 * L32 + y3 * L33
    return L

_y = L_function(_x, _x1, _x2, _x0, _x3, _y0, _y1, _y2, _y3)
print (_y)


plt.plot(_y)
plt.show()

Results:

a)

b)

c)

Splines

(Python code of the problem)

First of all, I used a 4x4 system equation solver. I obtained the values for x1, x2, x3 and x4. By using the formula, I evaluated the function and obtained the result of the function in point 27. The code of the last part can be seen below.

x1 = 14.263
x2 = -0.351354166667
x3 = 0.0054140625
x4 = -0.0000384114583333

x = 27

f27 = x1 + x2 * x + x3 * (x ** 2) + x4 * (x ** 3)

print f27

Results:

Integrals

(Python code of the problem using both methods)

from math import *
from pylab import *

def f(x):
    return 0.2 + 25 * x - 200 * (x **2) + 675 * (x ** 3) - 900 * (x ** 4) + 400 * (x ** 5)

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a + h
    for i in range(1,n/2 + 1):
        k += 4*f(x)
        x += 2*h

    x = a + 2*h
    for i in range(1,n/2):
        k += 2*f(x)
        x += 2*h
    return (h/3)*(f(a)+f(b)+k)

def trapezoidal(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

n = 2

def printFunction(x, n):
    print "   n       |     Value"
    for i in range(10):

        print "n = ", n, "      ", x
        n = n * 2

print " Simpson method: "
printFunction(simpson(f, 0.0, 0.8, n), n)

print " Trapezoidal method: "
printFunction(trapezoidal(f, 0.0, 0.8, n), n)

Conslcusion:

The trapezoidal rule is not as accurate as Simpson’s Rule when the underlying function is smooth, because Simpson’s rule uses quadratic approximations instead of linear approximations. The formula is usually given in the case of an odd number of equally spaced points.

Laboratory work No. 4

Matrices 1

(C++ code of the problem)

#include<iostream>

using namespace std;

int main(void){
    float a[10][10], b[10], x[10], y[10];
    int n = 3;
    int iterations = 20;

    //initial values
    x[0] = x[1] = x[2] = 0;
    //values for the matrix a
    a[0][0] = 1;
    a[0][1] = -2;
    a[1][0] = 2;
    a[1][1] = 1;


    //values from the rightpart of equal sign
    b[0] = -1;
    b[1] = 3;

    while (iterations > 0){
        for (int i = 0; i < n; i++){
            y[i] = (b[i] / a[i][i]);
            for (int j = 0; j < n; j++){
                if (j == i)
                    continue;
                y[i] = y[i] - ((a[i][j] / a[i][i]) * x[j]);
                x[i] = y[i];
            }
            cout << "x" << i+1 << " = " << y[i] << "\t\t\t";
        }
        cout << "\n";
        iterations--;
    }
    return 0;
}

Results:

x1 x2 x3
-1 5 -7.26705e-40
9 -15 .90044e-39
-31 65 -1.16082e-38
129 -255 4.64262e-38
-511 1025 -1.85711e-37
2049 -4095 7.42839e-37
-8191 16385 -2.97136e-36
32769 -65535 1.18854e-35
-131071 262145 -4.75418e-35
524289 -1.04858e+06 1.90167e-34
-2.09715e+06 4.1943e+06 -7.60668e-34
8.38861e+06 -1.67772e+07 3.04267e-33
-3.35544e+07 6.71089e+07 -1.21707e-32
1.34218e+08 -2.68435e+08 4.86828e-32
-5.36871e+08 1.07374e+09 -1.94731e-31
2.14748e+09 -4.29497e+09/td> 7.78924e-31
-8.58993e+09 1.71799e+10 -3.1157e-30
3.43597e+10 -6.87195e+10 1.24628e-2
-1.37439e+11 2.74878e+11 -4.98512e-29
5.49756e+11 -1.09951e+12 1.99405e-28

Conclusion:

As we analyze the results we observe that the method diverged. This reusults are due to the fact that the matrix does
not respect the conditions of the method.

Matrices 2

(C++ code of the problem)

#include<iostream>

using namespace std;

int main(void){
    float a[10][10], b[10], x[10], y[10];
    int n = 3;
    int iterations = 9;

    //initial values
    x[0] = x[1] = x[2] = 0;

    //values for the matrix a
    a[0][0] = 3;
    a[0][1] = -0.1;
    a[0][2] = -0.2;
    a[1][0] = 0.1;
    a[1][1] = 7;
    a[1][2] = -0.3;
    a[2][0] = 0.3;
    a[2][1] = -0.2;
    a[2][2] = 10;

    //values from the rightpart of equal sign
    b[0] = 7.85;
    b[1] = -19.3;
    b[2] = 71.4;

    while (iterations > 0){
        for (int i = 0; i < n; i++){
            y[i] = (b[i] / a[i][i]);
            for (int j = 0; j < n; j++){
                if (j == i)
                    continue;
                y[i] = y[i] - ((a[i][j] / a[i][i]) * x[j]);
                x[i] = y[i];
            }
            cout << "x" << i+1 << " = " << y[i] << "\t\t\t";
        }
        cout << "\n";
        iterations--;
    }
    return 0;
}

Results:

x1 = 0.0666667 x2 = 1.12048 x3 = -1.90959
x1 = -0.0899568 x2 = -4.1423 x3 = -2.01015
x1 = -0.272087 x2 = -4.41695 x3 = -2.01018
x1 = -0.281243 x2 = -4.4169 x3 = -2.0099
x1 = -0.281223 x2 = -4.41614 x3 = -2.00989
x1 = -0.281197 x2 = -4.4161 x3 = -2.00989
x1 = -0.281196 x2 = -4.4161 x3 = -2.00989
x1 = -0.281196 x2 = -4.4161 x3 = -2.00989
x1 = -0.281196 x2 = -4.4161 x3 = -2.00989
x1 = -0.281196 x2 = -4.4161 x3 = -2.00989

Conclusion:

Results:

x1 x2 x3
2.61667 -2.79452 7.00561
2.99056 -2.49962 7.00029
3.00003 -2.49999 7
3 -2.5 7
3 -2.5 7
3 -2.5 7
3 -2.5 7
3 -2.5 7
3 -2.5 7

Conclusion:

The second method converges, and it does it rapidly because the matrix satisfies the requirments.
As we see, a few iterations were sufficient to find out which is the solution.
Moreover, by analyzing information sources, we proved the method`s characteristics: Guass-Seidel method is one of the common methods employed for solving power flow equations.