python 3.x – I have written this code of a Finite Element Method for solving a Boundary value problem. But when I run this code it does not give me any output
Here is the code for the Boundary value problem and when I write this code in Matlab with the Matlab syntax it gives me a graph. But when I run this python code it does not gives me any result. Can anyone figure out that what is wrong with this code or what I am missing here.
import numpy as np
import matplotlib.pyplot as plt
def BVP2():
a=0 # interval [a,b]
b=1
n=10
h=(b-a)/n
nNode = n+1 # number of nodes
nElem = n # number of elements
x = np.zeros((nNode,1), dtype=int)
# DOmain discretization
for i in range(1, nNode):
x = a + (i-1)*h
print(x)
# Construction of finite elements
for i in range(1, nElem):
elem[i,:]= [i,i+1]
A = np.zeros((nNode,nNode))
B = np.zeros((nNode,nNode))
F = np.zeros((nNode,1))
u = np.zeros((nNode,1))
for i in range(1, nElem):
node1 = elem(i,1)
node2 = elem(i,2)
Ae = StiffMatrix1(x,h,node1,node2)
Be = StiffMatrix2(x,h,node1,node2)
Fe = Loadvector(x,h,node1,node2)
A[node1:node2,node1:node2] = A[node1:node2,node1:node2] + Ae
B[node1:node2,node1:node2] = B[node1:node2,node1:node2] + Be
F[node1:node2] = F[node1:node2] + Fe
# Implementation of BC's
Ared = A[2:nNode-1,2:nNode-1]
Bred = B[2:nNode-1,2:nNode-1]
Fred = F[2:nNode-1]
ured = u[2:nNode-1]
ured = [Ared + Bred]/Fred
u = np.block([[0],[ured], [0]])
# Exact Solution
uex = np.zeros((nNode,1))
for i in range(1, nNode):
uex = -(1/6)*x(i)**3 + (1/6)*x(i)
fig, ax = plt.subplots()
ax.plot(x, u)
plt.show()
def StiffMatrix1(x,h,node1,node2):
x1 = x(node1)
x2 = x(node2)
A11 = (1/h)
A12 = (-1/h)
A21 = (-1/h)
A22 = (1/h)
Ae = np.block([[np.c_[A11, A12]], [np.c_[A21, A22]]])
return Ae
def StiffMatrix2(x,h,node1,node2):
x1 = x(node1)
x2 = x(node2)
B11 = 0 #h/3
B12 = 0 #h/6
B21 = 0 #h/6
B22 = 0 #h/3
Be = np.block([[np.c_[B11, B12]], [np.c_[B21, B22]]])
return Be
def Loadvector(x,h,node1,node2):
x1 = x(node1)
x2 = x(node2)
Fe1 = h/6*(x2+2*x1)
Fe2 = h/6*(2*x2+x1)
Fe = np.block([[Fe1],[Fe2]])
return Fe
Your suggestions would be highly appreciated.
Read more here: Source link