Cómo podría convertir este código en python a matlab?
Publicado por Marisol (4 intervenciones) el 15/09/2020 20:22:02
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
%Este es el código en python, la verdad es que recién estoy aprendiendo un poco de matlab .Gracias de antemano.import numpy as np
import matplotlib.pyplot as pltimport pandas as pd
plt.style.use('classic')
import random
def sol_exacta(x,epsilon,b):
return (np.exp(b*x/epsilon)-1)/(np.exp(b/epsilon)-1)
def fem_ode(epsilon, beta, h, u_0, u_n):
n = int(1/h) + 1
Pe = beta*h/(2*epsilon) # Pe = k del problema
A = diagonalizar(-1-Pe, 2, Pe -1, n-2)
b = np.zeros(n-2)
b[ 0] = (1+Pe)*u_0
b[-1] = (1-Pe)*u_n
sol1 = metodo_thomas(A,b)
sol2 =np.append(u_0,sol1)
sol3=np.append(sol2,u_n)
return sol3
def diagonalizar(a,b,c,n):
D = [[a if i==j-1 else b if i==j else c if i==j+1 else 0 for i in range(n)] for j in range(n)]
return D
def metodo_thomas(A,b):
n = len(A)
if n == 1: return [b[0]/A[0][0]]
x = [0 for i in range(n)]
u = [A[i][i-1] for i in range(1,n)]
v = [A[i][i] for i in range(n)]
w = [A[i][i+1] for i in range(n-1)]
d = b.copy()
w[0] = w[0] / v[0]
d[0] = d[0] / v[0]
for i in range(1,n-1):
w[i] = w[i] / ( v[i] - w[i-1]*u[i])
for i in range(1,n):
d[i] = (d[i] - d[i-1]*u[i-1] ) / (v[i] - w[i-1]*u[i-1])
x[-1] = d[-1]
for i in reversed(range(n-1)):
x[i] = d[i] - w[i]*x[i+1]
return x
beta=float(input("beta : "));
h=float(input("paso : "))
u_0=float(input("condicion inicial : "))
u_n=float(input("condicion final : "))
n=int( input( "Numero de pruebas:" ))
epsilon = np.zeros(n+1, dtype=float)
Pe = np.zeros(n+1, dtype=float)
exac = np.zeros(n+1, dtype=float)
for i in range(1,n+1):
epsilon[i] = 10**(-i+int(n/2))
Pe[i] = beta*h/(2*epsilon[i])
print('iter: {} Cuando epsilon {} \t La constante de peckler esta dada por {} \t'.format(i,epsilon[i],Pe[i]) )
table1 = pd.DataFrame({'epsilon.': epsilon,'Constante Peckler.': Pe})
print(table1.to_latex(index=False)) #Esto nos servira para la documentacion
print('\n')
#print(table1)##DATOSm= int(1/h)+1
for i in range(1,n+1):
if (i!=n+1):
eps=epsilon[i]
print("Cuando epsilon es: {}".format(eps))
exac = sol_exacta(np.linspace(0,1,m),eps,1)
aprox = fem_ode(eps,beta,h,u_0,u_n)
table2 = pd.DataFrame({'Sol. Exacta.': exac,'Sol. Aprox.':aprox, ' ERROR ':np.abs(aprox-exac)})
print(table2.to_latex(index=False))
else:i+=1
number_of_colors = n+1
color = ["#"+''.join([random.choice('0123456789ABCDEF') for j in range(6)])
for i in range(number_of_colors)]
for i in range(1,n+1):
eps=epsilon[i]
print("Cuando epsilon es: {}".format(eps))
exac = sol_exacta(np.linspace(0,1,m),eps,1)
aprox = fem_ode(eps,beta,h,u_0,u_n)
plt.plot(exac,'r-',alpha=1,label='epsilon {}'.format(epsilon[i]), c=color[i] )
# plt.plot(aprox,'g^')plt.xlabel('tiempo')
plt.ylabel('solucion')
plt.title('Solucion Exacta')
plt.grid()
plt.legend(loc=0)
plt.savefig("solucion_exacta.jpg")
Valora esta pregunta
0
