Page MenuHomec4science
Files F97620366

darmadi.py
No OneTemporary
Actions

File Metadata

Created
Sun, Jan 5, 19:12

darmadi.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 10 12:16:05 2021
@author: ekinkubilay
"""
import numpy as np
import matplotlib.pyplot as plt
velocity = 1
t = 0.1
alpha = 1
upper = velocity**2*t/(4*alpha)
upper_limit = np.min([upper, 5])
lower_limit = 0
V = velocity/(2*alpha)
N = 20
points = mesh.coordinates()
source = np.zeros(np.shape(mesh.coordinates()))
source[:] = np.array([-1.8,0,4])
disp = source - points
R_squared = (disp*disp).sum(1)
R = np.sqrt(R_squared)
u = R*V
integral_value = np.zeros(len(points))
for i,j in enumerate(points):
integral_value[i] = integrate(function, u[i], 0.000001, upper_limit)[-1]
N = 200
def function(w, u):
return np.exp(-w-((u**2)/(4*w)))/(w**1.5)
def integrate(func, u, lower, upper):
integral = np.zeros(N-1)
delta = (upper-lower)/N
temp = 0
for i in range(N-1):
temp += (func(lower+(i+1)*delta, u)+func(lower+i*delta, u))*0.5*delta
integral[i] = temp
return integral
omega = np.linspace(0,5,N)
u = 1
u_range = np.array([3, 1, 0.5, 0.1])
for u in u_range:
f_omega = function(omega, u)
f_omega[0] = 0
#plt.plot(omega , f_omega, marker='.')
f_int = integrate(function, u, 0.00001, 5)
plt.plot(omega[:-1], f_int, marker='.', label='trapezoid')
print(f_int[-1])
simpson = integrate2(function, u, 0.00001, 5)
summation = np.zeros(N)
for i in range(N):
summation[i] = np.sum(simpson[0:i])
plt.plot(omega, summation, marker='x', label='simpson')
plt.legend()
def integrate2(func, u, lower, upper):
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
w = np.linspace(lower, upper, N)
integral_value = B.dot(func(u,w))
return integral_value*(upper-lower)/(3*N)
simpson = integrate2(function, 1, 0.000001, 5)
summation = np.zeros(N)
for i in range(N):
summation[i] = np.sum(simpson[0:i])
plt.plot(omega, summation)
trial = np.zeros(N-1)
for j,i in enumerate(omega[1:]):
trial[j] = integrate2(function, u, 0.000001, i)
plt.plot(omega[1:], trial)
N = 20
B = np.ones(N)
B[1:-1:2] = 4
B[2:-1:2] = 2
print(B)

Event Timeline

c4science · Help