热力学问题：运用牛顿拉弗逊迭代法确定相界
您已经了解了吉布斯自由能如何决定纯元素和混合物中的稳定相位。您还介绍了一些混合物如何理想，而另一些则不理想。非理想解的一个例子是正则解，其中正则解参数可能导致较低温度下的混性间隙。
您会记得，为了确定相图的稳定相位和相位边界，可以在多个相位可能稳定的情况下使用通用切线方法。从数学上讲，这可以通过分析方式对理想混合物进行，但实际上不可能以分析方式解决常规解决方案。在这种情况下，我们必须采用数值方法。
您将实现牛顿-拉夫森方法，以确定虚构的A-B合金系统中固体单相场和两相场之间的相边界，其中其中一个固相表现为正则解。
您应该提交两个文件：
单个pdf文档（*.pdf），包括
您所做的任何工作，都可以使用Netwon-Raphson方法将问题重新构成一个问题[30分]
显示500K至1500K之间的相位边界的相图，包括相关标签[20点]
简短的讨论（最多300字），讨论您的方法对起点的选择有多敏感，为什么会这样，以及如何使代码更健壮[20分]
您的Python实现（*.py），导致输出pdf文档中包含的图形。确保注释包含在您的代码[30分]中，并且代码可以执行。
这种虚观的A-B合金中每个纯相（以J/摩尔为单位）的吉布斯自由能可以根据以下方程建模：
图片
虽然α阶段表现得像一个理想的解决方案，β相位表现为正则解，具有正正则解参数：
图片2
提示：
您将看到这如何归结为求解一个由两个方程组来计算两个变量（相位边界的位置，这是B在α和β阶段）。我们简要介绍了如何处理多变量函数系统，但在这项任务中，您可以简单地组合两个方程来计算出β使用简单的牛顿-拉夫森方法的相位边界。然后，您可以使用另一个方程来计算出α相位边界。如果你愿意，你可以改用雅各比安来解决这个问题。两个阶段边界的位置将在两者之间Xb=0和Xb=1。如果您的解决方案在区间之外出现奇，Netwon-Raphson方法将无法收敛。如果您要解决的函数被写为Xb，那么Xb=0.999是牛顿-拉夫森方法的良好起。
各位能给点思路吗
https://blog.csdn.net/szlcw1/article/details/24327639/
需要公式的原理还是相关的程序
目前我将问题简化了，公式如下
1，为函数的原式，2为一阶导函数。
目前需要将原函数进行牛顿拉弗逊迭代，T为500-1500。需要根据最后迭代出来的根关于T（500-1500）之间进行作图。我在进行到牛顿拉弗逊迭代的时候出了点问题。显示如图
代码如下，我不知道是哪里出了问题，能有热心的大哥帮我分析分析吗，可以教教我之后的代码可以怎么写吗。求求了，救救孩子。
import math as exp
import math as ln
import numpy as np
T = 1500
Ω = 6350
R = 8.314
GBβ = -45000 - 0.18T -(3.4 10(-3)) T**2
GBα = -32000 - 0.375T -(7.7610*(-3)) *T2
GAα = -35000 - 0.25 T - 0.0201 T**2
GAβ = -38000 - 0.13 T - (510(-3)) *T2
def f(xBβ): #B in β
return GAα + RTln(1-exp((GBβ - GBα + Ω*(1-xBβ)2)/(RT) + ln(xBβ)) )- GAβ - RT*ln(1-xBβ) - Ω * xBβ2
def p(xBβ):
return -((RT)/(1-exp((GBβ - GBα + Ω(1-xBβ)2)/(RT) + ln(xBβ)))) * ((1/xBβ) + (Ω*(2xBβ - 2))/(RT)) * exp((GBβ - GBα + Ω(1-xBβ)2)/(RT) + ln(xBβ)) + (RT)/(1-xBβ) - 2 Ω xBβ
n=0
xBβ= 0.999
yold = xBβ - f(xBβ)/p(xBβ)
tol = 1e-6
error = 100
if f(xBβ) == 0:
print(xBβ)
else:
while error < tol:
xBβold = xBβ
xBβ = xBβ - f(xBβ)/p(xBβ)
if p(xBβ) == 0:
print('error')
break
n = n+1
xBβnew = xBβ
error = np.abs(xBβnew - xBβold)
if n> 100:
break
print(xBβ)
不行，看了半天好像理解的有点偏了
