Transformation

2010 days ago by JAVIER

import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D 
       
def T(u,w): return np.array([u^2-w^2,2*u*w]) 
       
def JT(u,w): return 4*(u*u + w*w) 
       
def fun(x,y): return y 
       
x = np.linspace(0,1,2) y = np.linspace(0,1,5) 
       
class Square: def __init__(self,p1,p2,id=0): self.id = id self.p1 = p1 self.p2 = p2 self.a = np.array(p2[0]-p1[0],float) self.b = np.array(p2[1]-p1[1],float) self.q1 = p1 self.q2 = p1+ np.array([self.a,0]) self.q3 = self.q2+ np.array([0,self.b]) self.q4 = self.q1+ np.array([0,self.b]) self.n = 10 self.dx = self.a/self.n self.dy = self.b/self.n self.P = np.vstack((self.q1,self.q2,self.q3,self.q4,self.q1)) self.center = (p1+p2)/2. self.Tcenter= T(self.center[0],self.center[1]) self.E = [] self.H = [] def generate_segments(self): for k in range(len(self.P)-1): x_vec = np.linspace(self.P[k,0], self.P[k+1,0], self.n) y_vec = np.linspace(self.P[k,1], self.P[k+1,1], self.n) points= np.array([x_vec,y_vec]) self.L.append(points) def generate_data(self): x0,y0 = self.p1 x1,y1 = self.p2 x_lin = np.linspace(x0,x1,self.n) y_lin = np.linspace(y0,y1,self.n) for i in range(self.n+1): y_vec = y0*np.ones(self.n) self.E.append(np.array([x_lin,y_vec])) Tx,Ty = T(x_lin,y_vec) self.H.append(np.array([Tx,Ty])) #-------------------------------------------- x_vec = x0*np.ones(self.n) self.E.append(np.array([x_vec,y_lin])) Tx,Ty = T(x_vec,y_lin) self.H.append(np.array([Tx,Ty])) #-------------------------------------------- y0 += self.dy x0 += self.dx def plot_original(self,axes1,col): for xvec,yvec in self.E: axes1.plot(xvec,yvec,color=col,linewidth=2) axes1.text(self.center[0], self.center[1], str(self.id), fontsize=15, fontweight='bold') def plot_transform(self,axes1,col): for xvec,yvec in self.H: axes1.plot(xvec,yvec,color=col,linewidth=2) axes1.text(self.Tcenter[0], self.Tcenter[1], str(self.id), fontsize=15, fontweight='bold') 
       
figure1 = plt.figure() axes1 = figure1.add_subplot(111) #--------------------------------- S = Square(np.array([1,1]),np.array([2,2])) S.generate_data() S.plot_original(axes1,[0,0,1]) plt.savefig('test.png') figure1.clear() 
       
figure2 = plt.figure() axes2 = figure2.add_subplot(111) #--------------------------------- S.plot_transform(axes2,[0,0,1]) plt.savefig('test2.png') figure1.clear() 
       
x0 = 1 y0 = 1 x1 = 2 y1 = 2 n = 4 x_mesh = np.linspace(x0,x1,n) y_mesh = np.linspace(y0,y1,n) L = [] A = 0 F = 0 count = -1 for i in range(n-1): y_0 = y_mesh[i] y_1 = y_mesh[i+1] for j in range(n-1): count += 1 x_0 = x_mesh[j] x_1 = x_mesh[j+1] du = x_1 - x_0 dw = y_1 - y_0 p1 = np.array([x_0,y_0]) p2 = np.array([x_1,y_1]) mid = (p1 + p2)/2 jac = JT(mid[0],mid[1]) JxW = du*dw*jac A += JxW rp = np.array([mid[0]^2 - mid[1]^2,2*mid[0]*mid[1]]) F += fun(rp[0],rp[1])*JxW print count, fun(rp[0],rp[1]) S = Square(p1,p2,count) S.generate_data() L.append(S) 
       
0 2.72222222222
1 3.5
2 4.27777777778
3 3.5
4 4.5
5 5.5
6 4.27777777778
7 5.5
8 6.72222222222
0 2.72222222222
1 3.5
2 4.27777777778
3 3.5
4 4.5
5 5.5
6 4.27777777778
7 5.5
8 6.72222222222
figure3 = plt.figure() axes3 = figure3.add_subplot(111) #--------------------------------- color_vector = [] for s in L: col = np.random.rand(3) color_vector.append(col) s.plot_original(axes3,col) plt.savefig('test3.png') figure3.clear() 
       
figure4 = plt.figure() axes4 = figure4.add_subplot(111) #--------------------------------- for k,s in enumerate(L): color_vector.append(col) s.plot_transform(axes4,color_vector[k]) plt.savefig('test4.png') figure4.clear() 
       
from scipy.integrate import dblquad as quad2d exact_area = quad2d(lambda u,w: 4*(u*u + w*w) , 1, 2, lambda w: 1, lambda w: 2) print 'exact area = ', exact_area print 'approx area =', A exact_integral = quad2d(lambda u,w: 2*u*w*4*(u*u + w*w) , 1, 2, lambda u: 1, lambda u: 2) print 'exact integral = ', exact_integral print 'approx integral =', F 
       
exact area =  (18.666666666666664, 2.0724163126336254e-13)
approx area = 18.5925925926
exact integral =  (90.0, 9.992007221626409e-13)
approx integral = 89.0
exact area =  (18.666666666666664, 2.0724163126336254e-13)
approx area = 18.5925925926
exact integral =  (90.0, 9.992007221626409e-13)
approx integral = 89.0
2*7/6*7/6 
       
49/18
49/18