非线性方程的数值解法----python

二分法

import math
def f(x):
	f = 2*x**3 - 5*x - 1
	return f

a = 1
b = 2
x0 = (a+b)/2
x1 , count= 0 ,0
while f(x0)!=0 and abs(x1-x0)>=1e-2 :
	if f(x0)*f(a)<0 :
		b = x0
	else:
		a = x0
	x1 = x0
	x0 = (a+b)/2
	count+=1
print(count)
print(x0)

迭代法

import math
def f(x):
	x = math.log(x+2,10)
	return x
x0, x1,count= 1, 0, 0
while abs(x1-x0)>=1e-6:
	x1 = f(x0)
	x0 = f(x1)
	count+=1
print(count)
print(x0)

Newton迭代法

import sympy as sy
import math
def f(x):
	y = x**3 + 2*x**2+10*x-20
	return y
x = sy.symbols('x')
y = x**3 + 2*x**2+10*x-20
z1 = sy.diff(y)
z2 = sy.diff(z1)

z3 = sy.diff(z2) #z3=6
a, b = 0, 2
#z1 = float(z1.evalf(subs={x:x}))
x0= -1
count = 0
x1 = 1.5
if float(y.evalf(subs={x:x1}))*float(z2.evalf(subs={x:x1}))>0:
	while abs(x1-x0)>=1e-6:
		x0 = x1
		t = float(y.evalf(subs={x:x0}))/float(z1.evalf(subs={x:x0}))
		print(t)
		x1 = x1 - t
		count = count + 1
		print(x1 ,'*')
print(count)
print(x1)

 

import numpy as np
import sympy as sy
import math
def f(x,p,q,m,n):
	f = x ** m - (1+ q/p)*x**(m-n) + q / p
	return f
x = sy.symbols('x')
p, q, m, n = 200, 2282, 639,420
y = x ** m - (1+ q/p)*x**(m-n) + q / p
z1 = sy.diff(y)
x0 = 1 
x1 = 1.5
while  abs(x1 - x0)>=1e-4:
	x0 = x1
	x1 = x1-float(y.evalf(subs={x:x0}))/float(z1.evalf(subs={x:x0}))
print(x1-1)

弦截法

def f(x):
	f = x ** 3 + 2 *x -6
	return f
x0 = 1.5
x1 = 2.0
count = 0
while abs(x1 - x0) > 1e-5:
	x2 = x1 - f(x1)*(x1 - x0)/(f(x1)-f(x0))
	x0 = x1
	x1 = x2
	count += 1
print(count)
print(x1)

scipy求解

import scipy.optimize as opt
import numpy as np
import math as m
def f1(x):
	return [m.exp(-m.exp(-(x[0]+x[1])))-x[1]- x[1]*(1 + x[0] ** 2),x[0] * m.cos(x[1]) + x[1] * m.sin(x[0])-1/2]

def jac(x):
	return np.array([[m.exp(-m.exp(-(x[0]+x[1])))*m.exp(-(x[0]+x[1]))-2*x[0]*x[1],m.cos(x[1])+x[1]*m.cos(x[1])],
		             [m.exp(-m.exp(-(x[0]+x[1])))*m.exp(-(x[0]+x[1]))-1-x[1]**2,-x[0]*m.sin(x[1])+m.sin(x[0])]])
sol = opt.root(f1,[0,0],jac = jac)
print(sol.x)
print(f1(sol.x))


import scipy.optimize as opt
import math
def f(x):
	f = 3*x**2 - math.exp(x)
	return f
print(opt.bisect(f, -1, 0))
print(opt.newton(f, 0, fprime=lambda x:6*x - math.exp(x)))