import pandas as pd #Importando librerias
import numpy as np
from math import *
import matplotlib.pyplot as plt
def biseccion(f,a,b,tol): #Comprobando si cumple el teorema de Bolzano
if f(a)*f(b) > 0:
print("Invalid interval does not satisfy Bolzano's theorem")
return
error = 1e3
X_anterior = 0
i = 1; I = []; Xa = []; Xb = []; Xm = []; #Creamos listas
Fa = []; Fb = []; E = [];Fx = []
while error > tol:
m = (a+b)/2
X_actual = m
error = (abs(X_actual-X_anterior)/X_actual) *100
I.append(i); Xa.append(a); Xb.append(b);
Xm.append(m); Fa.append(f(a)); Fb.append(f(b)); #Rellenamos las listas
E.append(error);Fx.append(f(m))
if f(a)*f(m) < 0:
b = m
else:
a = m
X_anterior = X_actual
i += 1
d = {"I":I,"F(x)":Fx,"Xa":Xa,"Xb":Xb,"Xm":Xm,"Fa":Fa,"Fb":Fb,"E %":E,}
TT = pd.DataFrame(d)
TT.set_index("I",inplace=True)
print(TT)
f = lambda x: x**2 - 2
biseccion(f,1,2,0.01)
F(x) Xa Xb Xm Fa Fb E %
I
1 0.250000 1.000000 2.000000 1.500000 -1.000000 2.000000 100.000000
2 -0.437500 1.000000 1.500000 1.250000 -1.000000 0.250000 20.000000
3 -0.109375 1.250000 1.500000 1.375000 -0.437500 0.250000 9.090909
4 0.066406 1.375000 1.500000 1.437500 -0.109375 0.250000 4.347826
5 -0.022461 1.375000 1.437500 1.406250 -0.109375 0.066406 2.222222
6 0.021729 1.406250 1.437500 1.421875 -0.022461 0.066406 1.098901
7 -0.000427 1.406250 1.421875 1.414062 -0.022461 0.021729 0.552486
8 0.010635 1.414062 1.421875 1.417969 -0.000427 0.021729 0.275482
9 0.005100 1.414062 1.417969 1.416016 -0.000427 0.010635 0.137931
10 0.002336 1.414062 1.416016 1.415039 -0.000427 0.005100 0.069013
11 0.000954 1.414062 1.415039 1.414551 -0.000427 0.002336 0.034518
12 0.000263 1.414062 1.414551 1.414307 -0.000427 0.000954 0.017262
13 -0.000082 1.414062 1.414307 1.414185 -0.000427 0.000263 0.008632
def plotf(f,a,b,n):
x = np.linspace(a,b, num=n)
y = f(x)
fig, ax = plt.subplots()
ax.plot(x,y)
ax.grid()
ax.axhline(y=0, color='r')
ax.axvline(x=0, color='r')
plotf(f,0,2,100)
def falsapos(f,a,b,tol):
if f(a)*f(b) > 0:
print("Invalid interval does not satisfy Bolzano's theorem")
return
error = 1e3
X_anterior = 0;
i = 1; I = []; Xa = []; Xb = []; Xm = [];
Fa = []; Fb = []; Fm = []; E = [];
while error > tol:
m = a - f(a)*(b-a)/(f(b)-f(a))
X_actual = m
error = (abs(X_actual - X_anterior)/X_actual)*100
I.append(i); Xa.append(a); Xb.append(b);
Xm.append(m); Fa.append(f(a)); Fb.append(f(b));
Fm.append(f(m)); E.append(error);
if f(a)*f(m) < 0:
b = m
else:
a = m
X_anterior = X_actual
i += 1
d = {"I":I,"Xa":Xa,"Xb":Xb,"Xm":Xm,"Fa":Fa,"Fm":Fm,"Fb":Fb,"E":E}
TT = pd.DataFrame(d)
TT.set_index("I",inplace=True)
print(TT.to_string())
def puntof(f,x0,tol):
xold = x0
e = 100
a = lambda x: x
i = 1
I = []
Xr = []
Fx = []
Gx = []
Ea = []
Xr.append(x0)
Fx.append(f(x0))
Ea.append(0)
while e > tol:
xr = f(xold) + a(xold)
e = (abs(xr - xold)/xr) * 100
Gx.append(xr)
xold = xr
I.append(i)
Xr.append(xr)
Fx.append(f(xr))
Ea.append(e)
i += 1
Gx.append(f(xold) + a(xold))
d = {"Xr":Xr,"G(x)":Gx,"F(x)":Fx,"E":Ea}
TT = pd.DataFrame(d)
print(TT.to_string())
f1 = lambda x: exp(-x) - x
x0 = 0
tol = 0.1
puntof(f1,x0,tol)
Xr G(x) F(x) E
0 0.000000 1.000000 1.000000 0.000000
1 1.000000 0.367879 -0.632121 100.000000
2 0.367879 0.692201 0.324321 171.828183
3 0.692201 0.500474 -0.191727 46.853639
4 0.500474 0.606244 0.105770 38.309147
5 0.606244 0.545396 -0.060848 17.446790
6 0.545396 0.579612 0.034217 11.156623
7 0.579612 0.560115 -0.019497 5.903351
8 0.560115 0.571143 0.011028 3.480867
9 0.571143 0.564879 -0.006264 1.930804
10 0.564879 0.568429 0.003549 1.108868
11 0.568429 0.566415 -0.002014 0.624419
12 0.566415 0.567557 0.001142 0.355568
13 0.567557 0.566909 -0.000648 0.201197
14 0.566909 0.567276 0.000367 0.114256
15 0.567276 0.567068 -0.000208 0.064752
def derivada(f,x0):
h=0.000000001
d=(f(x0+h)-f(x0))/h
if abs(d)<0.000001:
d=0
return(d)
def newton(f,x0,tol):
e=100
n=1;i=[];X0=[];X=[];h=[];E=[];
if derivada(f,x0)==0:
print("Enter another closest value")
else:
while e>tol:
x=x0-(f(x0))/(derivada(f,x0))
e=(abs(x-x0)/x) * 100
fa=f(x0)
i.append(n);X0.append(x0);X.append(x);
h.append(fa); E.append(e);
x0=x
n+=1
d = {"i":i,"x0":X0,"x":X,"f(x)":h,"E %":E}
TT = pd.DataFrame(d)
TT.set_index("i",inplace=True)
print(TT.to_string())
f2 = lambda x : exp(-x) - x
x0 = 0
newton(f2,x0,0.01)
x0 x f(x) E %
i
1 0.000000 0.500000 1.000000e+00 100.000000
2 0.500000 0.566311 1.065306e-01 11.709289
3 0.566311 0.567143 1.304504e-03 0.146728
4 0.567143 0.567143 1.964546e-07 0.000022
# Importing NumPy Library
import numpy as np
import sys
# Reading number of unknowns
n = int(input('Enter number of unknowns: '))
# Making numpy array of n x n+1 size and initializing
# to zero for storing augmented matrix
a = np.zeros((n,n+1))
# Making numpy array of n size and initializing
# to zero for storing solution vector
x = np.zeros(n)
# Reading augmented matrix coefficients
print('Enter Augmented Matrix Coefficients:')
for i in range(n):
for j in range(n+1):
a[i][j] = int(input( 'a['+str(i)+']['+ str(j)+']='))
# Applying Gauss Elimination
for i in range(n):
if a[i][i] == 0.0:
sys.exit('Divide by zero detected!')
for j in range(i+1, n):
ratio = a[j][i]/a[i][i]
for k in range(n+1):
a[j][k] = a[j][k] - ratio * a[i][k]
# Back Substitution
x[n-1] = a[n-1][n]/a[n-1][n-1]
for i in range(n-2,-1,-1):
x[i] = a[i][n]
for j in range(i+1,n):
x[i] = x[i] - a[i][j]*x[j]
x[i] = x[i]/a[i][i]
# Displaying solution
print('\nRequired solution is: ')
for i in range(n):
print('X%d = %0.2f' %(i,x[i]), end = '\t')
Enter Augmented Matrix Coefficients:
Required solution is:
X0 = 1.00 X1 = 3.00 import numpy as np
def diy_lu(a):
"""Construct the LU decomposition of the input matrix.
Naive LU decomposition: work column by column, accumulate elementary triangular matrices.
No pivoting.
"""
N = a.shape[0] # gives number of filas
u = a.copy() # Copy of a
L = np.eye(N) # gives a matriz of N with 1 in his diagonal
for j in range(N-1):
lam = np.eye(N) # Matriz of 1 in diagonal
gamma = u[j+1:, j] / u[j, j] # gamma = [u(1,0)/u(0,0)]
lam[j+1:, j] = -gamma
u = lam @ u
print("U", j,"\n", u,"\n")
lam[j+1:, j] = gamma
L = L @ lam
return L, u
a = np.array([[10,2,-1],
[-3,-6,2],
[1,1,5]])
l, u = diy_lu(a)
print("este es L","\n",l,"\n")
print("este es U","\n",u,"\n")
U 0
[[10. 2. -1. ]
[ 0. -5.4 1.7]
[ 0. 0.8 5.1]]
U 1
[[10. 2. -1. ]
[ 0. -5.4 1.7 ]
[ 0. 0. 5.35185185]]
este es L
[[ 1. 0. 0. ]
[-0.3 1. 0. ]
[ 0.1 -0.14814815 1. ]]
este es U
[[10. 2. -1. ]
[ 0. -5.4 1.7 ]
[ 0. 0. 5.35185185]]
# Gauss Seidel Iteration
# Defining equations to be solved
# in diagonally dominant form
f1 = lambda x,y,z: (7.85 + 0.1*y+0.2*z)/3
f2 = lambda x,y,z: (-19.3-0.1*x+0.3*z)/7
f3 = lambda x,y,z: (71.4-0.3*x+0.2*y)/10
# Initial setup
x0 = 0
y0 = 0
z0 = 0
count = 1
# Reading tolerable error
e = float(input('Enter tolerable error: '))
# Implementation of Gauss Seidel Iteration
print('\nCount\tx\ty\tz\n')
condition = True
while condition:
x1 = f1(x0,y0,z0)
y1 = f2(x1,y0,z0)
z1 = f3(x1,y1,z0)
print('%d\t%0.4f\t%0.4f\t%0.4f\n' %(count, x1,y1,z1))
e1 = abs(x0-x1);
e2 = abs(y0-y1);
e3 = abs(z0-z1);
count += 1
x0 = x1
y0 = y1
z0 = z1
condition = e1>e and e2>e and e3>e
print('\nSolution: x=%0.3f, y=%0.3f and z = %0.3f\n'% (x1,y1,z1))
Count x y z
1 2.6167 -2.7945 7.0056
2 2.9906 -2.4996 7.0003
Solution: x=2.991, y=-2.500 and z = 7.000