Page MenuHomec4science
Files F88829061

untitled1.py
No OneTemporary
Actions

File Metadata

Created
Sun, Oct 20, 21:36

untitled1.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 3 22:01:47 2022
@author: kubilay
"""
import numpy as np
import matplotlib.pyplot as plt
N = 20
def function(w, u):
return 2*np.exp(-w-((u**2)/(4*w)))/(w**1.5)
def integrate2(func, u, lower, upper):
B = np.ones(400)
B[1:-1:2] = 4
B[2:-1:2] = 2
limit = 1e-1
w = np.linspace(lower, limit, 400)
integral_value = B.dot(func(u,w))
w = np.linspace(lower, upper, 400)
integral_value = B.dot(func(u,w))
return integral_value*(upper-lower)/(3*400)
simpson = integrate2(function, 1, 0.000001, 5)
int_value = np.zeros(N)
func_value = np.zeros(N)
int_value_adaptive = np.zeros(N)
omega = np.linspace(0.000001, 1000, N)
for i,f in enumerate(omega):
func_value[i] = function(f, 0.10)
int_value[i] = integrate2(function, 0.1, 0.000001, f)
int_value_adaptive[i] = integrate2(function, 0.1, 0.000001, f/2) + integrate2(function, 0.1, f/2, f)
plt.plot(omega, func_value, marker=".")
plt.show()
plt.plot(omega,int_value)
plt.plot(omega, int_value_adaptive, marker=".")
plt.plot(omega, int_value - int_value_adaptive, marker=".")
plt.show()
def integrate_all(function, u, lower, upper):
omega = np.linspace(lower, upper, 40)
int_value = np.zeros(len(omega))
for i,f in enumerate(omega):
int_value[i] = integrate2(function, u, min(omega), f)
return int_value
def integrate_all_adaptive(function, u, omega):
int_value = np.zeros(len(omega))
quad_1 = np.zeros(len(omega))
quad_2 = np.zeros(len(omega))
int_value = integrate_all(function, u, min(omega), max(omega))
quad_1 = integrate_all(function, u, min(omega), 0.2*max(omega))
quad_2 = integrate_all(function, u, 0.2*max(omega), max(omega))
while np.any(np.abs(int_value-quad_1-quad_2)>7):
int_value = quad_1 + quad_2
quad_1 = integrate_all_adaptive(function, u, np.linspace(min(omega),max(omega)*0.2, 40))
#quad_2 = integrate_all_adaptive(function, u, np.linspace(max(omega/2),max(omega), 40))
#plt.plot(int_value-quad_1-quad_2)
return quad_1 + quad_2
omega = np.linspace(0.0000001, 1000, 40)
u = 0.1
normal_integration = integrate_all(function, u, min(omega), max(omega))
adaptive_integration = integrate_all_adaptive(function, u, omega)
plt.plot(omega, normal_integration)
plt.plot(omega, adaptive_integration)
#%%
u = 1e-5
v = 0.8
N = 1000
alpha = 1
def new_function(u, v, t):
return np.exp(-(u**2)/(4*alpha*t) - (t*v**2)/(4*alpha))/(t**1.5)
value = np.zeros(N)
for i,t in enumerate(np.linspace(1e-10,1,N)):
value[i] = new_function(u,v,t)
plt.plot(value, marker=".")
u /= np.sqrt(4*alpha)
v /= np.sqrt(4*alpha)
t_max = (np.sqrt(16*u**2*v**2 + 9)+3)/(4*v**2)
from scipy.optimize import fsolve
func = lambda t : (u**2/(1))*(1/t - 1/t_max)-((v**2)/(1))*(t_max-t) - np.log((0.001)*(t_max/t)**1.5)
plt.plot(np.linspace(0,100,N), func(np.linspace(0,100,N)))
plt.xlim([-1,1])
plt.ylim([-1,1])
plt.plot(np.linspace(0,100,N), value, marker=".")
t_initial_guess = t_max - 1
t_solution = fsolve(func, t_initial_guess)
t_solution
t_initial_guess = t_max + 10
t_solution = fsolve(func, t_initial_guess)
t_solution
#%%
u = 1e-3
v = 8
N = 10000
alpha = 1e-2
t_max = (np.sqrt(u**2*v**2/(alpha**2) + 9)+3)/((v**2)/alpha)
def f_1(u,v,t,alpha):
return -u**2/(4*alpha*t) - t*v**2/(4*alpha)
def f_2(u,v,t,alpha):
t_max = (np.sqrt(u**2*v**2/(alpha**2) + 9)+3)/(alpha*v**2)
return f_1(u,v,t_max, alpha) + np.log(0.01*(t/t_max)**1.5)
def f_3(t,u,v,alpha):
u = args[0]
v = args[1]
alpha = args[2]
return f_2(u,v,t,alpha)-f_1(u,v,t,alpha)
t = np.linspace(1e-10, 40, N)
plt.plot(t, f_1(u,v,t,alpha), marker=".")
plt.plot(t, f_2(u,v,t,alpha), marker=".")
plt.plot(t, f_2(u,v,t,alpha)-f_1(u,v,t,alpha), marker=".")
plt.xlim([25,30])
plt.ylim([-1,1])
plt.show()
from scipy.optimize import newton
u_range = np.linspace(1e-10, 0.1, N)
t_upper_limit = []
t_lower_limit = []
for u in u_range:
args=(u,v,alpha)
t_upper_limit.append(newton(f_3, 1, args=(u,v,alpha)))
#t_lower_limit.append(fsolve(f_3, 0, args=(u,v,alpha)))
plt.plot(u_range, t_upper_limit, marker=".")
plt.plot(u_range, t_lower_limit)
plt.xlim([0,0.0001])
delta_square = np.einsum('ij,ij->i', mesh.coordinates()-scan.initial_loc-scan.velocity*(t-scan.t_start), mesh.coordinates()-scan.initial_loc-scan.velocity*(t-scan.t_start))
delta = np.sqrt(delta_square)
A_min, B_min, A_max, B_max = get_interpolation_limits(0.00377/2, v, alpha)
dof = len(delta_square)
N = 10
time_ranges = np.zeros((dof, N))
time_range = np.linspace(0, 1, N)
for i in range(dof):
time_ranges[i] = np.linspace(A_min*delta[i]+B_min, A_max*delta[i]+B_max, N)
lower = 0.00001
upper = t
A_min, B_min, A_max, B_max = get_interpolation_limits(0.00377/4, velocity_mag, alpha)
delta = np.sqrt(delta_square)
N = 100
dof = len(delta_square)
time_ranges = np.zeros((dof, N))
for i in range(dof):
time_ranges[i] = np.linspace(max([A_min*delta[i]+B_min, lower]), min([A_max*delta[i]+B_max, upper]), N)
time_range = np.linspace(lower, upper, N)

Event Timeline

c4science · Help