科学计算课程题目

在一维中求解非线性方程

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)

这是我根据老师的课程写出来的课后题，如果存在错误，请大家指出，十分感谢！！！