Home

Education & employment

Publications

Presentations

Research

CV

 

A program for diffusion-couple fitting

Equation for fitting:

c=c(x): the experimentally measured concentration profile.

A:concentration in material M.

C: concentration in material N.

H: interface position.

B=4Dt: D is the diffusion coefficient and t is the annealing duration.

Program:

 

Code:

'====================================================================
'程序名称:多元非线性拟合-高斯-牛顿迭代法
'程序版本:1.0.0.1
'编程语言:vb
'编程环境:Microsoft Visual Basic 6.0
'作 者: 费宏展 Fei Hongzhan
'完成日期:2011.03.05
'====================================================================
Const maxdatanumber = 500
Const p_max = 4 。
Const Pi = 3.14159
Const p = 4

Dim caled As Boolean
Dim n As Integer
Dim X(0 To maxdatanumber) As Double
Dim Y(0 To maxdatanumber) As Double
Dim Rr(0 To p_max) As Double。
Dim x_down_temp As Integer
Dim y_down_temp As Integer
Dim x_up_temp As Integer
Dim y_up_temp As Integer
Dim xmax, xmin, ymax, ymin As Double
Dim xmax_ori, xmin_ori, ymax_ori, ymin_ori As Double
Dim XaxisWidth, YaxisHeight As Integer
'-------------------------------------------------------------------
Private Const OFN_EXPLORER = &H80000
Private Const OFN_ALLOWMULTISELECT = &H200
Private Declare Function GetOpenFileName Lib "comdlg32.dll" Alias "GetOpenFileNameA" (pOpenfilename As OPENFILENAME) As Long
Private Type OPENFILENAME
lStructSize As Long
hWndOwner As Long
hInstance As Long
lpstrFilter As String
lpstrCustomFilter As String
nMaxCustFilter As Long
nFilterIndex As Long
lpstrFile As String
nMaxFile As Long
lpstrFileTitle As String
nMaxFileTitle As Long
lpstrInitialDir As String
lpstrTitle As String
flags As Long
nFileOffset As Integer
nFileExtension As Integer
lpstrDefExt As String
lCustData As Long
lpfnHook As Long
lpTemplateName As String
End Type
'------------------------------------------------------------------
'====================================================================
Private Sub Button_readdata_Click()
Me.Cls
caled = False
n = 0
Dim arr
Dim str As String
On Error GoTo fileopenerror
Open filename_text.Text For Input As #1
Do While Not EOF(1)
Line Input #1, str
n = n + 1
arr = Split(str, vbTab)
X(n) = arr(0)
Y(n) = arr(1)
Loop
Close #1
data_number.Text = "共有" & n & "行数据。"
Call paintfig(X(), Y())
Exit Sub
fileopenerror:
Close #1
n = 0
MsgBox ("文件不存在或者格式不符合要求!")
Exit Sub
End Sub
'====================================================================
Function findmax(ByRef p() As Double) As Double
Dim i As Integer
findmax = p(1)
For i = 1 To n
If findmax < p(i) Then
findmax = p(i)
End If
Next i
End Function
'====================================================================
Function findmin(ByRef p() As Double) As Double
Dim i As Integer
findmin = p(1)
For i = 1 To n
If findmin > p(i) Then
findmin = p(i)
End If
Next i
End Function
'====================================================================
Private Function paintfig(ByRef X() As Double, ByRef Y() As Double)
Dim deltaX, deltay As Double
Dim i As Integer
XaxisWidth = 8100 - 600
YaxisHeight = 7000 - 2800
xmax = findmax(X())
xmin = findmin(X())
ymax = findmax(Y())
ymin = findmin(Y())
deltaX = (xmax - xmin)
deltay = (ymax - ymin)
xmax = xmax + 0.05 * deltaX
xmin = xmin - 0.05 * deltaX
ymax = ymax + 0.05 * deltay
ymin = ymin - 0.05 * deltay
Text_xmin.Text = xmin
Text_xmax.Text = xmax
Text_ymin.Text = ymin
Text_ymax.Text = ymax
xmin_ori = xmin
xmax_ori = xmax
ymin_ori = ymin
ymax_ori = ymax
ForeColor = vbBlue
DrawWidth = 5
For i = 1 To n - 1
Me.PSet (600 + (X(i) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (Y(i) - ymin) * (YaxisHeight) / (ymax - ymin))), vbBlue
'Line (600 + (x(i) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp1 - ymin) * (YaxisHeight) / (ymax - ymin)))-(600 + (x(i + 1) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp2 - ymin) * (YaxisHeight) / (ymax - ymin)))
Next
End Function
'====================================================================
Private Sub Cal_Click()
Call Disableall

Dim YY(0 To maxdatanumber) As Double
Dim DD(0 To maxdatanumber, 0 To p_max) As Double
Dim DDXX(0 To p_max, 0 To p_max) As Double
Dim DDXX_YY(0 To p_max) As Double
Dim DDXXop(0 To p_max, 0 To p_max) As Double
Dim Rr0(0 To p_max) As Double
Dim caltime As Integer
Dim SSR, tempSSR As Double
Dim ERR As Double
Dim i, j, k As Integer
Dim a(0 To p_max + 2, 0 To 2 * p + 3) As Double
Dim m(0 To 2 * p_max + 3) As Double
Dim s As Integer
'result_text.Text = "迭代次数" & vbTab & "A" & vbTab & "B" & vbTab & "C" & vbTab & "D" & vbTab & "误差率" & vbTab & "相对误差" & vbCrLf
result_text.Text = "迭代次数 A B C D 误差率 相对误差" & vbCrLf
On Error GoTo errormsg
Rr(1) = A_start.Text
Rr(2) = B_start.Text
Rr(3) = C_start.Text
Rr(4) = D_start.Text

caltime = cal_time_number.Text
ERR = 1
SSR = 0
For s = 1 To caltime '迭代的循环
tempSSR = SSR '记录上一次迭代的SSR
SSR = 0 '
For i = 1 To n
YY(i) = Y(i) - resulttemp(Rr(1), Rr(2), Rr(3), Rr(4), X(i)) '得到YY矩阵
SSR = SSR + Abs(YY(i)) '初始值计算出来的SSR
DD(i, 1) = (1 - erf((X(i) - Rr(4)) / Sqr(Rr(2)))) / 2
DD(i, 2) = (Rr(1) - Rr(3)) * (X(i) - Rr(4)) * Exp(-(X(i) - Rr(4)) * (X(i) - Rr(4)) / Rr(2)) / Rr(2) / Sqr(Pi * Rr(2)) / 2
DD(i, 3) = (1 + erf((X(i) - Rr(4)) / Sqr(Rr(2)))) / 2
DD(i, 4) = (Rr(1) - Rr(3)) * Exp(-(X(i) - Rr(4)) * (X(i) - Rr(4)) / Rr(2)) / Sqr(Pi * Rr(2))
'对系数求偏导
Next
'-----------------------------------------------------------
For i = 0 To n
DD(i, 0) = 0
Next i
'计算DDXX,矩阵DD与DD的转置矩阵相乘

For i = 0 To p
For j = 0 To p
DDXX(i, j) = 0
Next j
Next i

For i = 1 To p
For j = 1 To p
For k = 1 To n
DDXX(i, j) = DDXX(i, j) + (DD(k, i) * DD(k, j))
Next k
Next j
Next i

DDXX(0, 0) = n
For i = 1 To p
For k = 1 To n
DDXX(0, i) = DDXX(0, i) + DD(k, i)
DDXX(i, 0) = DDXX(i, 0) + DD(k, i)
Next k
Next i

'计算DDXY_YY矩阵相乘
For i = 0 To p
DDXX_YY(i) = 0
Next i
For i = 1 To n
DDXX_YY(0) = DDXX_YY(0) + YY(i)
Next i
For i = 1 To p
For j = 1 To n
DDXX_YY(i) = DDXX_YY(i) + YY(j) * DD(j, i)
Next j
Next i

'计算DDXXop-逆矩阵

For i = 1 To p + 1
For j = 1 To p + 1
a(i, j) = DDXX(i - 1, j - 1)
Next j
Next i

For i = 1 To p + 1
For j = p + 2 To 2 * p + 2
If j - i = p + 1 Then
a(i, j) = 1
Else
a(i, j) = 0
End If
Next j
Next i

For k = 1 To p
For i = k + 1 To p + 1
m(i) = a(i, k) / a(k, k)
For j = k To 2 * p + 2
a(i, j) = a(i, j) - m(i) * a(k, j)
Next j
Next i
Next k

For k = p + 1 To 2 Step -1
For i = k - 1 To 1 Step -1
For j = p + 2 To 2 * p + 2
a(i, j) = a(i, j) - a(i, k) * a(k, j) / a(k, k)
Next j
Next i
Next k

For i = 1 To p + 1
For j = i + 1 To 2 * p + 2
a(i, j) = a(i, j) / a(i, i)
Next j
Next i

For i = 0 To p
For j = 0 To p
DDXXop(i, j) = a(i + 1, j + p + 2)
Next j
Next i

'计算系数修正值Rr0()矩阵。
For i = 0 To p
Rr0(i) = 0
Next i

For i = 0 To p
For j = 0 To p
Rr0(i) = Rr0(i) + DDXXop(i, j) * DDXX_YY(j)
Next j

Rr(i) = Rr(i) + Rr0(i)
Next i
'--------------------------------------------------
If fix_A.Value = 1 Then
Rr(1) = A_start.Text
End If
If fix_B.Value = 1 Then
Rr(2) = B_start.Text
End If
If fix_C.Value = 1 Then
Rr(3) = C_start.Text
End If
If fix_D.Value = 1 Then
Rr(4) = D_start.Text
End If

SSR = SSR / n '计算平均偏差值
ERR = Abs(tempSSR - SSR) / SSR '计算误差率。

Dim Rrtemp(0 To p_max) As Double
Dim SSRtemp, ERRtemp As Double
For i = 1 To p
Rrtemp(i) = Format(Rr(i), "#0.0000000000")
Next
SSRtemp = Format(SSR, "#0.00000000")
ERRtemp = Format(ERR, "#0.00000000")

result_text.Text = result_text.Text & s & vbTab & Left(str(Rrtemp(1)), 8) & vbTab & Left(str(Rrtemp(2)), 9) & vbTab & Left(str(Rrtemp(3)), 8) & vbTab & Left(str(Rrtemp(4)), 7) & vbTab & Left(str(ERRtemp * 100), 5) & "%" & vbTab & Left(str(SSRtemp), 7) & vbCrLf
Next
ForeColor = vbRed
DrawWidth = 3
Dim temp1, temp2 As Double
For i = 1 To n - 1
temp1 = (Rr(3) - Rr(1)) * erf((X(i) - Rr(4)) / Sqr(Rr(2))) / 2 + (Rr(1) + Rr(3)) / 2
temp2 = (Rr(3) - Rr(1)) * erf((X(i + 1) - Rr(4)) / Sqr(Rr(2))) / 2 + (Rr(1) + Rr(3)) / 2
Line (600 + (X(i) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp1 - ymin) * (YaxisHeight) / (ymax - ymin)))-(600 + (X(i + 1) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp2 - ymin) * (YaxisHeight) / (ymax - ymin)))
Next

Call Enableall

caled = True

Exit Sub

errormsg:
MsgBox "初始值或迭代次数选择不合理"
Call Enableall
End Sub
'====================================================================
Function resulttemp(AA As Double, BB As Double, CC As Double, DD As Double, xxii As Double) As Double '从初始值计算结果
resulttemp = (AA + CC) / 2 - (AA - CC) / 2 * erf((xxii - DD) / Sqr(BB))
End Function
'====================================================================
Private Sub Disableall()
file_open.Enabled = False
Button_readdata.Enabled = False
startdata.Enabled = False
cal_time.Enabled = False
filename_text.Enabled = False
data_number.Enabled = False
Result_as_start_A.Enabled = False
Result_as_start_B.Enabled = False
Result_as_start_C.Enabled = False
Result_as_start_D.Enabled = False
A_start.Enabled = False
B_start.Enabled = False
C_start.Enabled = False
D_start.Enabled = False
cal_time_number.Enabled = False
Cal.Enabled = False
clear_data.Enabled = False
result_text.Enabled = False
fig_return.Enabled = False
End Sub
'====================================================================
Private Sub Enableall()
file_open.Enabled = True
Button_readdata.Enabled = True
startdata.Enabled = True
cal_time.Enabled = True
filename_text.Enabled = True
data_number.Enabled = True
Result_as_start_A.Enabled = True
Result_as_start_B.Enabled = True
Result_as_start_C.Enabled = True
Result_as_start_D.Enabled = True
A_start.Enabled = True
B_start.Enabled = True
C_start.Enabled = True
D_start.Enabled = True
cal_time_number.Enabled = True
Cal.Enabled = True
clear_data.Enabled = True
result_text.Enabled = True
fig_return.Enabled = True
End Sub
'====================================================================
Private Sub clear_data_Click()
result_text.Text = ""
Me.Cls
End Sub
'====================================================================
Private Sub Label_email_Click()
MsgBox ("feihongzhan@gmail.com")
End Sub
'====================================================================
Private Sub file_open_Click()
Debug.Print FileOpen
End Sub
'====================================================================
Function FileOpen() As String
Dim OpenFN As OPENFILENAME
Dim Rc As Long
With OpenFN
.hWndOwner = 0
.hInstance = App.hInstance
.lpstrTitle = "" 'Title
.lpstrFilter = "*.*" 'Filter
.nFilterIndex = 1 'FilterIndex
.lpstrInitialDir = "" 'StartDir
.lpstrFile = FileOpen & String$(255, Chr$(0))
.nMaxFile = 255
.lpstrFileTitle = .lpstrFile
.nMaxFileTitle = 255
.flags = OFN_ALLOWMULTISELECT Or OFN_EXPLORER
.lStructSize = Len(OpenFN)
End With
Rc = GetOpenFileName(OpenFN)
If Rc Then
FileOpen = Left$(OpenFN.lpstrFile, OpenFN.nMaxFile)
filename_text.Text = FileOpen
End If
End Function
'====================================================================
Function erf(ByRef X As Double) As Double
Dim k, i As Integer
Dim deltaX As Double
k = 100
deltaX = X / k
erf = 0
For i = 0 To k - 1
erf = erf + deltaX * Exp((-1) * (i * deltaX + 0.5 * deltaX) * (i * deltaX + 0.5 * deltaX))
Next
erf = erf * (2 / Sqr(Pi))
End Function
'====================================================================
Private Sub Result_as_start_A_Click()
A_start.Text = Rr(1)
End Sub
Private Sub Result_as_start_B_Click()
B_start.Text = Rr(2)
End Sub
Private Sub Result_as_start_C_Click()
C_start.Text = Rr(3)
End Sub
Private Sub Result_as_start_D_Click()
D_start.Text = Rr(4)
End Sub
'====================================================================
Private Sub startdata_Click()
A_start.Text = Rr(1)
B_start.Text = Rr(2)
C_start.Text = Rr(3)
D_start.Text = Rr(4)
End Sub
'====================================================================
Private Sub Text_ymax_KeyPress(KeyAscii As Integer)
If KeyAscii = 13 Then
ymax = Text_ymax.Text
Call paintagain
End If
End Sub
Private Sub Text_ymin_KeyPress(KeyAscii As Integer)
If KeyAscii = 13 Then
ymin = Text_ymin.Text
Call paintagain
End If
End Sub
Private Sub Text_xmax_KeyPress(KeyAscii As Integer)
If KeyAscii = 13 Then
xmax = Text_xmax.Text
Call paintagain
End If
End Sub
Private Sub Text_xmin_KeyPress(KeyAscii As Integer)
If KeyAscii = 13 Then
xmin = Text_xmin.Text
Call paintagain
End If
End Sub
Private Sub Form_MouseDown(Button As Integer, Shift As Integer, xxx As Single, yyy As Single)
If Button = 1 Then
x_down_temp = xxx
y_down_temp = yyy
End If
End Sub
Private Sub Form_Mouseup(Button As Integer, Shift As Integer, xxx As Single, yyy As Single)
Dim tempex As Integer
Dim x_up_temp, y_up_temp As Integer
If Button = 1 Then
x_up_temp = xxx
y_up_temp = yyy
End If
If x_down_temp >= 600 And x_down_temp <= 8100 And x_up_temp >= 600 And x_up_temp <= 8100 And y_down_temp >= 2800 And y_down_temp <= 7000 And y_up_temp >= 2800 And y_up_temp <= 7000 Then
If x_up_temp < x_down_temp Then
tempex = x_up_temp
x_up_temp = x_down_temp
x_down_temp = tempex
End If
If y_up_temp > y_down_temp Then
tempex = y_up_temp
y_up_temp = y_down_temp
y_down_temp = tempex
End If
If x_up_temp = x_down_temp Then
x_up_temp = 8100
x_down_temp = 600
End If
If y_up_temp = y_down_temp Then
y_up_temp = 2800
y_down_temp = 7000
End If
ymax = ymax - (Text_ymax.Text - Text_ymin.Text) * (y_up_temp - 2800) / (7000 - 2800)
ymin = ymax - (Text_ymax.Text - Text_ymin.Text) * (y_down_temp - y_up_temp) / (7000 - 2800)
xmax = xmin + (Text_xmax.Text - Text_xmin.Text) * (x_up_temp - 600) / (8100 - 600)
xmin = xmax - (Text_xmax.Text - Text_xmin.Text) * (x_up_temp - x_down_temp) / (8100 - 600)
Call paintagain
Text_xmax.Text = xmax
Text_ymax.Text = ymax
Text_xmin.Text = xmin
Text_ymin.Text = ymin
End If
End Sub
Private Sub fig_return_Click()
xmin = xmin_ori
xmax = xmax_ori
ymin = ymin_ori
ymax = ymax_ori
Text_xmin.Text = xmin
Text_xmax.Text = xmax
Text_ymin.Text = ymin
Text_ymax.Text = ymax
Call paintagain
End Sub
'====================================================================
Private Sub paintagain()
On Error GoTo errormsg
Me.Cls
Dim i As Integer
ForeColor = vbBlue
DrawWidth = 5
For i = 1 To n - 1
Me.PSet (600 + (X(i) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (Y(i) - ymin) * (YaxisHeight) / (ymax - ymin))), vbBlue
Next
If caled = True Then
ForeColor = vbRed
DrawWidth = 3
Dim temp1, temp2 As Double
For i = 1 To n - 1
temp1 = (Rr(3) - Rr(1)) * erf((X(i) - Rr(4)) / Sqr(Rr(2))) / 2 + (Rr(1) + Rr(3)) / 2
temp2 = (Rr(3) - Rr(1)) * erf((X(i + 1) - Rr(4)) / Sqr(Rr(2))) / 2 + (Rr(1) + Rr(3)) / 2
Line (600 + (X(i) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp1 - ymin) * (YaxisHeight) / (ymax - ymin)))-(600 + (X(i + 1) - xmin) * (XaxisWidth) / (xmax - xmin), (7000 - (temp2 - ymin) * (YaxisHeight) / (ymax - ymin)))
Next
End If
Exit Sub
errormsg:
End Sub
'====================================================================

   
>