In [1]:

    

import numpy as np
from matplotlib import pylab as plt
%matplotlib inline




num_of_points = 201
h = np.pi / (num_of_points - 1)
num_of_elements = num_of_points - 1
N = 3 * num_of_elements  # global size of matrixes
x = np.linspace(0, np.pi, num_of_points)
X = np.linspace(0, np.pi, N)

exact_solution = np.zeros(num_of_points)
K_j = np.zeros((4, 4))  # local
M_j = np.zeros((4, 4))  # local
K = np.zeros([N, N])  # global
M = np.zeros([N, N])  # global
matr = np.zeros([N, N])  # K+M
I = np.zeros(N)


# exact solution
def solution(x):
    return 1 / 110 * (
            (np.exp(-np.sqrt(10) * x) * (np.sqrt(10) * np.exp(np.sqrt(10) * np.pi) + np.sqrt(10) * np.exp(
                2 * np.sqrt(10) * x) - 10 * np.exp(np.sqrt(10) * x) * np.sin(x) + 10 * np.exp(
                np.sqrt(10) * np.pi + np.sqrt(10) * x) * np.sin(x))) /
            (np.exp(np.sqrt(10) * np.pi) - 1))


# values of the integrals
def f1(x1, x2):
    c = (x2 - x1) / 2
    b = (x2 + x1) / 2
    return (np.cos(b) * ((27 - 13 * c * c) * np.sin(c) + c * (4 * c * c - 27) * np.cos(c)) + c * np.sin(b) * (
            (4 * c * c - 9) * np.sin(c) + 9 * c * np.cos(c))) / (4 * c * c * c * c)


def f2(x1, x2):
    c = (x2 - x1) / 2
    b = (x2 + x1) / 2
    return 9 * (3 * np.cos(b) * ((c * c - 3) * np.sin(c) + 3 * c * np.cos(c)) + c * np.sin(b) * (
            np.sin(c) - c * np.cos(c))) / (4 * c * c * c * c)


def f3(x1, x2):
    c = (x2 - x1) / 2
    b = (x2 + x1) / 2
    return -9 * (3 * np.cos(b) * ((c * c - 3) * np.sin(c) + 3 * c * np.cos(c)) + c * np.sin(b) * (
            -np.sin(c) + c * np.cos(c))) / (4 * c * c * c * c)


def f4(x1, x2):
    c = (x2 - x1) / 2
    b = (x2 + x1) / 2
    return (np.cos(b) * ((27 * c - 4 * c * c * c) * np.cos(c) + (13 * c * c - 27) * np.sin(c)) + c * np.sin(b) * (
            (4 * c * c - 9) * np.sin(c) + 9 * c * np.cos(c))) / (4 * c * c * c * c)


    

In [2]:

    

for i in np.arange(num_of_points):
    exact_solution[i] = solution(x[i])
#matrix filling
K_j[0][0] = K_j[3][3] = 1.85
K_j[0][1] = K_j[1][0] = K_j[3][2] = K_j[2][3] = -2.3625
K_j[1][1] = K_j[2][2] = 5.4
K_j[0][2] = K_j[2][0] = K_j[3][1] = K_j[1][3] = 0.675
K_j[0][3] = K_j[3][0] = -0.1625
K_j[1][2] = K_j[2][1] = -3.7125

M_j[0][0] = M_j[3][3] = 16 / 105
M_j[0][1] = M_j[1][0] = M_j[3][2] = M_j[2][3] = 33 / 280
M_j[1][1] = M_j[2][2] = 27 / 35
M_j[0][2] = M_j[2][0] = M_j[3][1] = M_j[1][3] = -3 / 70
M_j[0][3] = M_j[3][0] = 19 / 840
M_j[1][2] = M_j[2][1] = -27 / 280

for i in np.arange(4):
    for k in np.arange(4):
        K_j[i][k] = K_j[i][k] * 2 / h
        M_j[i][k] = M_j[i][k] * 5 * h
        K[i][k] = K_j[i][k]
        M[i][k] = M_j[i][k]

for s in range(3 * num_of_elements - 3):
    if (s % 3 == 0):
        M[2 + s:6 + s, 2 + s:6 + s] += M_j
        K[2 + s:6 + s, 2 + s:6 + s] += K_j
ind = 0
for i in range(num_of_elements - 1):
    I[ind] = f2(x[i], x[i + 1])
    I[ind + 1] = f3(x[i], x[i + 1])
    I[ind + 2] = f4(x[i], x[i + 1]) + f1(x[i], x[i + 1])
    ind += 3

for i in range(N):
    I[i] = I[i] * h / 2

for i in range(N):
    for k in range(N):
        matr[i][k] = K[i][k] + M[i][k]
vect = np.linalg.solve(matr, I)
#print(I)
#print(M_j)


    

    




[1.85056948e-05 7.40203245e-05 3.08414991e-05 1.11027148e-04
 1.66528080e-04 9.25168876e-05 2.03521207e-04 2.58994748e-04
 1.54169449e-04 2.95965050e-04 3.51397512e-04 2.15783971e-04
 3.88335868e-04 4.43713575e-04 2.77345252e-04 4.80610870e-04
 5.35920157e-04 3.38838102e-04 5.72767289e-04 6.27994510e-04
 4.00247349e-04 6.64782386e-04 7.19913914e-04 4.61557841e-04
 7.56633458e-04 8.11655690e-04 5.22754450e-04 8.48297841e-04
 9.03197202e-04 5.83822078e-04 9.39752920e-04 9.94515864e-04
 6.44745656e-04 1.03097613e-03 1.08558914e-03 7.05510152e-04
 1.12194496e-03 1.17639457e-03 7.66100575e-04 1.21263697e-03
 1.26690974e-03 8.26501974e-04 1.30302978e-03 1.35711232e-03
 8.86699446e-04 1.39310108e-03 1.44698005e-03 9.46678138e-04
 1.48282866e-03 1.53649076e-03 1.00642325e-03 1.57219037e-03
 1.62562237e-03 1.06592004e-03 1.66116417e-03 1.71435287e-03
 1.12515384e-03 1.74972810e-03 1.80266039e-03 1.18411002e-03
 1.83786031e-03 1.89052312e-03 1.24277403e-03 1.92553906e-03
 1.97791940e-03 1.30113142e-03 2.01274271e-03 2.06482766e-03
 1.35916776e-03 2.09944974e-03 2.15122645e-03 1.41686876e-03
 2.18563877e-03 2.23709445e-03 1.47422016e-03 2.27128852e-03
 2.32241049e-03 1.53120782e-03 2.35637787e-03 2.40715351e-03
 1.58781767e-03 2.44088581e-03 2.49130260e-03 1.64403576e-03
 2.52479150e-03 2.57483700e-03 1.69984820e-03 2.60807424e-03
 2.65773609e-03 1.75524124e-03 2.69071348e-03 2.73997943e-03
 1.81020119e-03 2.77268882e-03 2.82154672e-03 1.86471451e-03
 2.85398004e-03 2.90241783e-03 1.91876773e-03 2.93456709e-03
 2.98257282e-03 1.97234753e-03 3.01443008e-03 3.06199190e-03
 2.02544066e-03 3.09354930e-03 3.14065548e-03 2.07803408e-03
 3.17190523e-03 3.21854415e-03 2.13011475e-03 3.24947854e-03
 3.29563869e-03 2.18166986e-03 3.32625010e-03 3.37192009e-03
 2.23268667e-03 3.40220095e-03 3.44736951e-03 2.28315258e-03
 3.47731235e-03 3.52196835e-03 2.33305519e-03 3.55156579e-03
 3.59569819e-03 2.38238214e-03 3.62494292e-03 3.66854085e-03
 2.43112127e-03 3.69742566e-03 3.74047835e-03 2.47926055e-03
 3.76899611e-03 3.81149294e-03 2.52678812e-03 3.83963662e-03
 3.88156711e-03 2.57369223e-03 3.90932976e-03 3.95068355e-03
 2.61996134e-03 3.97805833e-03 4.01882523e-03 2.66558400e-03
 4.04580537e-03 4.08597532e-03 2.71054897e-03 4.11255417e-03
 4.15211725e-03 2.75484516e-03 4.17828826e-03 4.21723471e-03
 2.79846161e-03 4.24299141e-03 4.28131164e-03 2.84138761e-03
 4.30664768e-03 4.34433221e-03 2.88361253e-03 4.36924134e-03
 4.40628089e-03 2.92512596e-03 4.43075696e-03 4.46714238e-03
 2.96591766e-03 4.49117935e-03 4.52690167e-03 3.00597756e-03
 4.55049362e-03 4.58554401e-03 3.04529577e-03 4.60868511e-03
 4.64305494e-03 3.08386262e-03 4.66573949e-03 4.69942027e-03
 3.12166857e-03 4.72164266e-03 4.75462609e-03 3.15870429e-03
 4.77638084e-03 4.80865877e-03 3.19496065e-03 4.82994051e-03
 4.86150499e-03 3.23042869e-03 4.88230847e-03 4.91315170e-03
 3.26509968e-03 4.93347180e-03 4.96358617e-03 3.29896505e-03
 4.98341786e-03 5.01279594e-03 3.33201645e-03 5.03213434e-03
 5.06076889e-03 3.36424574e-03 5.07960922e-03 5.10749316e-03
 3.39564493e-03 5.12583078e-03 5.15295724e-03 3.42620631e-03
 5.17078761e-03 5.19714990e-03 3.45592232e-03 5.21446864e-03
 5.24006025e-03 3.48478563e-03 5.25686307e-03 5.28167768e-03
 3.51278913e-03 5.29796044e-03 5.32199195e-03 3.53992589e-03
 5.33775063e-03 5.36099309e-03 3.56618924e-03 5.37622381e-03
 5.39867148e-03 3.59157268e-03 5.41337048e-03 5.43501784e-03
 3.61606995e-03 5.44918148e-03 5.47002318e-03 3.63967502e-03
 5.48364798e-03 5.50367888e-03 3.66238204e-03 5.51676148e-03
 5.53597663e-03 3.68418543e-03 5.54851379e-03 5.56690846e-03
 3.70507980e-03 5.57889709e-03 5.59646673e-03 3.72505999e-03
 5.60790388e-03 5.62464417e-03 3.74412109e-03 5.63552701e-03
 5.65143380e-03 3.76225838e-03 5.66175965e-03 5.67682903e-03
 3.77946739e-03 5.68659534e-03 5.70082359e-03 3.79574387e-03
 5.71002795e-03 5.72341155e-03 3.81108381e-03 5.73205169e-03
 5.74458735e-03 3.82548342e-03 5.75266113e-03 5.76434576e-03
 3.83893915e-03 5.77185120e-03 5.78268190e-03 3.85144768e-03
 5.78961714e-03 5.79959125e-03 3.86300592e-03 5.80595458e-03
 5.81506964e-03 3.87361102e-03 5.82085949e-03 5.82911324e-03
 3.88326036e-03 5.83432819e-03 5.84171860e-03 3.89195158e-03
 5.84635736e-03 5.85288261e-03 3.89968251e-03 5.85694403e-03
 5.86260250e-03 3.90645125e-03 5.86608558e-03 5.87087588e-03
 3.91225613e-03 5.87377976e-03 5.87770071e-03 3.91709572e-03
 5.88002468e-03 5.88307531e-03 3.92096883e-03 5.88481879e-03
 5.88699835e-03 3.92387450e-03 5.88816090e-03 5.88946886e-03
 3.92581201e-03 5.89005020e-03 5.89048622e-03 3.92678088e-03
 5.89048622e-03 5.89005020e-03 3.92678088e-03 5.88946886e-03
 5.88816090e-03 3.92581201e-03 5.88699835e-03 5.88481879e-03
 3.92387450e-03 5.88307531e-03 5.88002468e-03 3.92096883e-03
 5.87770071e-03 5.87377976e-03 3.91709572e-03 5.87087588e-03
 5.86608558e-03 3.91225613e-03 5.86260250e-03 5.85694403e-03
 3.90645125e-03 5.85288261e-03 5.84635736e-03 3.89968251e-03
 5.84171860e-03 5.83432819e-03 3.89195158e-03 5.82911324e-03
 5.82085949e-03 3.88326037e-03 5.81506964e-03 5.80595458e-03
 3.87361102e-03 5.79959125e-03 5.78961714e-03 3.86300592e-03
 5.78268190e-03 5.77185120e-03 3.85144768e-03 5.76434576e-03
 5.75266113e-03 3.83893915e-03 5.74458735e-03 5.73205169e-03
 3.82548342e-03 5.72341155e-03 5.71002795e-03 3.81108381e-03
 5.70082359e-03 5.68659534e-03 3.79574387e-03 5.67682903e-03
 5.66175965e-03 3.77946739e-03 5.65143380e-03 5.63552701e-03
 3.76225839e-03 5.62464417e-03 5.60790388e-03 3.74412109e-03
 5.59646673e-03 5.57889709e-03 3.72506000e-03 5.56690846e-03
 5.54851379e-03 3.70507980e-03 5.53597663e-03 5.51676148e-03
 3.68418544e-03 5.50367888e-03 5.48364798e-03 3.66238205e-03
 5.47002318e-03 5.44918148e-03 3.63967502e-03 5.43501784e-03
 5.41337048e-03 3.61606996e-03 5.39867148e-03 5.37622381e-03
 3.59157268e-03 5.36099309e-03 5.33775063e-03 3.56618925e-03
 5.32199195e-03 5.29796044e-03 3.53992591e-03 5.28167768e-03
 5.25686307e-03 3.51278914e-03 5.24006025e-03 5.21446864e-03
 3.48478564e-03 5.19714990e-03 5.17078761e-03 3.45592233e-03
 5.15295724e-03 5.12583078e-03 3.42620632e-03 5.10749316e-03
 5.07960922e-03 3.39564494e-03 5.06076889e-03 5.03213434e-03
 3.36424574e-03 5.01279594e-03 4.98341786e-03 3.33201646e-03
 4.96358617e-03 4.93347180e-03 3.29896506e-03 4.91315170e-03
 4.88230847e-03 3.26509969e-03 4.86150499e-03 4.82994051e-03
 3.23042870e-03 4.80865877e-03 4.77638084e-03 3.19496065e-03
 4.75462609e-03 4.72164266e-03 3.15870429e-03 4.69942027e-03
 4.66573949e-03 3.12166858e-03 4.64305494e-03 4.60868511e-03
 3.08386263e-03 4.58554401e-03 4.55049362e-03 3.04529579e-03
 4.52690167e-03 4.49117935e-03 3.00597757e-03 4.46714238e-03
 4.43075696e-03 2.96591766e-03 4.40628089e-03 4.36924134e-03
 2.92512597e-03 4.34433221e-03 4.30664768e-03 2.88361254e-03
 4.28131164e-03 4.24299141e-03 2.84138762e-03 4.21723471e-03
 4.17828826e-03 2.79846163e-03 4.15211725e-03 4.11255417e-03
 2.75484516e-03 4.08597532e-03 4.04580537e-03 2.71054898e-03
 4.01882523e-03 3.97805833e-03 2.66558401e-03 3.95068355e-03
 3.90932976e-03 2.61996135e-03 3.88156711e-03 3.83963662e-03
 2.57369225e-03 3.81149294e-03 3.76899611e-03 2.52678812e-03
 3.74047835e-03 3.69742566e-03 2.47926056e-03 3.66854085e-03
 3.62494292e-03 2.43112128e-03 3.59569819e-03 3.55156579e-03
 2.38238215e-03 3.52196835e-03 3.47731235e-03 2.33305520e-03
 3.44736951e-03 3.40220095e-03 2.28315259e-03 3.37192009e-03
 3.32625010e-03 2.23268668e-03 3.29563869e-03 3.24947854e-03
 2.18166987e-03 3.21854415e-03 3.17190523e-03 2.13011477e-03
 3.14065548e-03 3.09354930e-03 2.07803408e-03 3.06199190e-03
 3.01443008e-03 2.02544069e-03 2.98257282e-03 2.93456709e-03
 1.97234754e-03 2.90241783e-03 2.85398004e-03 1.91876774e-03
 2.82154672e-03 2.77268882e-03 1.86471452e-03 2.73997943e-03
 2.69071348e-03 1.81020119e-03 2.65773609e-03 2.60807424e-03
 1.75524125e-03 2.57483700e-03 2.52479150e-03 1.69984822e-03
 2.49130260e-03 2.44088581e-03 1.64403577e-03 2.40715351e-03
 2.35637787e-03 1.58781769e-03 2.32241049e-03 2.27128852e-03
 1.53120782e-03 2.23709445e-03 2.18563877e-03 1.47422017e-03
 2.15122645e-03 2.09944974e-03 1.41686877e-03 2.06482766e-03
 2.01274271e-03 1.35916778e-03 1.97791940e-03 1.92553906e-03
 1.30113143e-03 1.89052312e-03 1.83786031e-03 1.24277403e-03
 1.80266039e-03 1.74972810e-03 1.18411003e-03 1.71435287e-03
 1.66116417e-03 1.12515385e-03 1.62562237e-03 1.57219037e-03
 1.06592006e-03 1.53649076e-03 1.48282866e-03 1.00642327e-03
 1.44698005e-03 1.39310108e-03 9.46678138e-04 1.35711232e-03
 1.30302978e-03 8.86699460e-04 1.26690974e-03 1.21263697e-03
 8.26501988e-04 1.17639457e-03 1.12194496e-03 7.66100589e-04
 1.08558914e-03 1.03097613e-03 7.05510167e-04 9.94515864e-04
 9.39752920e-04 6.44745656e-04 9.03197202e-04 8.48297841e-04
 5.83822092e-04 8.11655690e-04 7.56633458e-04 5.22754464e-04
 7.19913914e-04 6.64782386e-04 4.61557855e-04 6.27994510e-04
 5.72767289e-04 4.00247363e-04 5.35920157e-04 4.80610870e-04
 3.38838102e-04 4.43713575e-04 3.88335868e-04 2.77345266e-04
 3.51397512e-04 2.95965050e-04 2.15783986e-04 2.58994748e-04
 2.03521207e-04 1.54169463e-04 1.66528080e-04 1.11027148e-04
 9.25169019e-05 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[[ 0.01196797  0.00925648 -0.00336599  0.0017765 ]
 [ 0.00925648  0.06058786 -0.00757348 -0.00336599]
 [-0.00336599 -0.00757348  0.06058786  0.00925648]
 [ 0.0017765  -0.00336599  0.00925648  0.01196797]]

In [4]:

    

plt.subplot(221)
plt.plot(X, vect, label="our_solution")
plt.legend()
plt.show()


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: