In [1]:
from IPython.display import display
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
k1 = 10
k2 = 1
a = 0
b = 0.5
c = 1
ua = 0
uc = 1
f = 1
In [3]:
h = 0.1
n = int((c - a) / h)
m = int((b-a)/c * n)-1
In [4]:
def mat_f(i, j):
if i == m:
if i == j:
return k1 + k2
if j == i-1:
return -k1
if j == i+1:
return -k2
if i < m:
cons = k1
elif i > m:
cons = k2
if i == j:
return 2*cons
elif j == i-1 or j == i+1:
return -cons
else:
return 0
In [5]:
M = np.fromfunction(np.vectorize(mat_f), (n-1, n-1))
M *= 1/h
In [6]:
y = h * f * np.ones(n-1)
y[0] += k1 * 1/h * ua
y[n-2] += k2 * 1/h * uc
In [7]:
x = np.dot(np.linalg.inv(M),y)
u = np.hstack((np.array([ua]), x, np.array([uc])))
In [8]:
plt.plot(u)
plt.show()
In [9]:
koef1 = k1 * np.ones(m+1)
koef2 = k2 * np.ones(n-m-1)
koeff = np.hstack([koef1, koef2])
udot = np.diff(u)/h
plt.plot(udot*koeff)
plt.show()
In [ ]: