在一维中求解非线性方程
1.课件01课后计算题
(1)
Solve the equation f
(
x
) =
x
3
−
2
x
2
= 0 from the initial guess
x0
= 1
.
3 and stopping tolerance
tol
= 10
−
6
. How many iterations
does Newton take to converge?
方程为f(x) = x^3 - 2x^2,初始猜测为x0=1.3,停止迭代的容许误差tol = 1e-6
下边为使用colab进行编程的代码
def myNewton(fnon, dfdx, x0, tol, maxk, *fnonargs):
# fnon - name of the nonlinear function f(x)
# dfdx - name of the derivative function df(x)
# x0 - initial guess for the solution x0
# tol - stopping tolerance for Newton's iteration
# maxk - maximum number of Newton iterations before stopping
# fnonargs - optional arguments that will be passed to the nonlinear function (useful for additional function parameters)
k = 0
x = x0
f = eval(fnon)(x,*fnonargs)
print(' k xk f(xk)')
# Main Newton loop
while (abs(f) > tol and k < maxk):
# Evaluate function derivative
d = eval(dfdx)(x,*fnonargs)
# Take Newton step
x = x - f/d
f = eval(fnon)(x,*fnonargs)
k += 1
print('{0:2.0f} {1:2.8f} {2:2.2e}'.format(k, x, abs(f)))
if (k == maxk):
print('Not converged')
else:
print('Converged')
#牛顿迭代法的方程
def squareRoot(x):
# Define function f(x) = x^3 - 2 * x^2
return (x*x*x - 2*x*x)
#原方程
def dSquareRoot(x):
# Define function df(x) = 3*x^2 - 4*x
return (3*x*x - 4*x)
#原方程的导数
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
#绘制部分的代码
x=np.arange(1, 2.25, 0.005)
#x的范围和范围的步数
y=squareRoot(x)
#y为原本的方程
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
#绘制图片的范围
ax.plot(x,y)
ax.plot(x,0*y)
ax.set_title('Function $f(x) = x^3 - 2x^2$')
#绘制出图的表头
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
myNewton('squareRoot', 'dSquareRoot', 1.3, 1e-6, 100)
#绘制结果(原始函数,原始函数的导数,初始值,停止迭代的容许误差,迭代次数)
下边是运算产生的图片和结果
图片:
结果:
(2)
Solve the equation f
(
x
) =
x
4
−
10
x
2
+ 0
.
1 = 0 from difffferent initial
guesses x
0
∈
[2
,
3]. How many difffferent real solutions can you find
that the method converges to?
方程为f(x) = x^4 - 10x^2 + 0.1,初始猜测的取值范围为[2,3],有多少可能得收敛步骤数
下边为colab的编程代码
#第一步的牛顿迭代函数代码与第一题相同
def squareRoot(x):
# Define function f(x) = x^4 - 10*x^2 + 0.1
return (x*x*x*x - 10*x*x + 0.1)
#原函数
def dSquareRoot(x):
# Define function df(x) = 4*x^2 - 20*x
return (4*x*x - 20*x)
#原函数求导
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
x=np.arange(2, 4, 0.005)
#x的取值范围,根据画出的图后期修改的范围
y=squareRoot(x)
fig=plt.figure()
ax=fig.add_axes([0,0,1,1])
ax.plot(x,y)
ax.plot(x,0*y)
ax.set_title('Function $f(x) = x^4 - 10*x^2 + 0.1$')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
#绘制图像的函数
myNewton('squareRoot', 'dSquareRoot', 3, 1e-6, 100)
# myNewton('squareRoot', 'dSquareRoot', 2.5, 1e-6, 100)
# myNewton('squareRoot', 'dSquareRoot', 2, 1e-6, 100)
#三种不同的初始值取值
下边为图片与结果
图片
结果
(3)
Consider the solution of the differential equation x
0
(
t
) =
−
x
(
t
)
2
from the initial value x
(0) = 1. We discretise this equation using
finite differences to obtain the discrete nonlinear equation:
F(
x
k
) :=
x
k
−
1
−
x
k
+ ∆
t
(
−
x
k
2
) = 0
Use Newton’s method to take 10 time steps with a constant
time-step ∆t
= 0
.
1. What is the numerical value of
x
(1)?
方程为F(x)= x_prev - x + dt * (-x**2),初始值为x0=1
下边是编译出的python代码
def F(x_prev, x, dt):
return x_prev - x + dt * (-x**2)
#原方程(Xk-1 - Xk + dt * (-Xk**2) )
def dF_dx(x, dt):
return -1 - 2 * dt * x
#原方程的导数
def newton_method(x0, dt, max_iter=4000, tol=1e-6):
x_prev = x0
for k in range(max_iter):
F_value = F(x_prev, x_prev, dt)
dF_dx_value = dF_dx(x_prev, dt)
if abs(dF_dx_value) < tol:
print(f"Derivative is close to zero at iteration {k}. Newton's method may not converge.")
raise ValueError("Derivative is close to zero. Newton's method may not converge.")
x = x_prev - F_value / dF_dx_value
if abs(x - x_prev) < tol:
return x
x_prev = x
print(f"Newton's method did not converge within {max_iter} iterations.")
raise ValueError("Newton's method did not converge within the specified iterations.")
#牛顿迭代法
def solve_differential_equation(initial_value, dt, num_steps):
x = initial_value
for k in range(num_steps):
x = newton_method(x, dt)
return x
#在设置的步骤中进行迭代
# 设置初始条件和时间步长
initial_value = 1.0
dt = 0.1
num_steps = 10
# 使用牛顿法解决微分方程
result = solve_differential_equation(initial_value, dt, num_steps)
print(f"The numerical value of x(1) is: {result}")
下边是编译出来的结果
The numerical value of x(1) is: 0.003152394061705063
牛顿迭代法中我认为比较重要的部分是
Xk+1 = Xk - f(Xk) / f'(Xk)
这是我根据老师的课程写出来的课后题,如果存在错误,请大家指出,十分感谢!!!