非线性方程？

[ForMyLove]

怎么用程序实现解非线性方程（组）？

[block]

要看你是要解决什么方面的问题，在非线性问题中同样的算法不同的问题得到的效果是不同的，如果精度要求不高使用对分法即可，工程实际应用可以采用拟牛顿算法，我在一些应用中（主要是多刚体运动分析中）使用拟牛顿算法求解过超过一千多个约束方程的非线性方程组，收敛效果非常好。

[block]

to:Montaque(Rainman)
冒昧的问一下，高斯消元算法处理线性问题有很多实例，但是它可以处理非线性问题吗？

[ForMyLove]

∑yi-n*a-b∑exp(-r*xi)=0
∑yi*exp(-r*xi)-a*∑exp(-r*xi)-b*∑exp(-2*r*xi)=0
b*∑xi*yi*exp(-r*xi)-b*a∑exp(-r*xi)-b*b∑xi*EXP(-2*r*xi)=0

[ForMyLove]

Q(a,b,r)=∑(yi-a-b*Exp(-r*xi))^2
上式中的yi和xi是从数据库中读取的，现在要求的是a,b,r 三个未知数。

[block]

Q(a,b,r)=∑(yi-a-b*Exp(-r*xi))^2中既然有三个未知数，那么应该有三个独立方程啊，否则没有定解！

[block]

根据你上面的要求应该给出类似下面这样的方程组，然后用拟牛顿算法即可处理。
f(1)=Q1(a,b,r)-∑(yi1-a-b*Exp(-r*xi1))^2
f(2)=Q2(a,b,r)-∑(yi2-a-b*Exp(-r*xi2))^2
f(3)=Q3(a,b,r)-∑(yi3-a-b*Exp(-r*xi3))^2

f(1)-f(3)，Q1(a,b,r)-Q1(a,b,r)，yi1-yi3，xi1-xi3均为已知量.

[block]

下面的内容给你留言时候提示留言太长，只好发到这里了。
下面的是实际使用的拟牛顿法解非线性方程组代码，你应该去掉里面捕捉错误用的自定义例外处理代码，如何你的程序中使用固定的方程组过程，可以不用定义委托，直接把委托的过程放到拟牛顿法解非线性方程组的过程中调用即可。
Public Delegate Sub pnlfsDelegate(ByRef x() As Double, ByRef y() As Double, ByRef xio() As Double)
'定义委托
Public Sub pnlfs(ByVal N As Integer, ByRef x() As Double, ByRef y() As Double, ByVal subguoceng As pnlfsDelegate,

ByVal xio() As Double)
'拟牛顿法解非线性方程组
Dim A(N, N) As Double, p As Double, h As Double, emax As Double
Dim w As Double, eps As Double
Dim i As Integer, i0 As Integer, i1 As Integer, i2 As Integer, j As Integer, k As Integer
Dim yesno As Boolean = True

h = 0.5
w = 0.5
eps = Val(1 / (10 ^ rlmaina.jsjdszBar.Value))
emax = 1.0E+300
Try
Call subguoceng(x, y, xio)

Do
Call subguoceng(x, y, xio)

yesno = True
For i = 1 To N
If System.Math.Abs(y(i)) > eps Then
yesno = False
For i0 = 1 To N
If System.Math.Abs(y(i0)) > emax Then
pgkeyexception.Source = "<错误讯息>：error:gyzcx.pnlfs" + Chr(10)
Throw pgkeyexception
Exit Sub
End If
Next i0
For j = 1 To N
x(j) = x(j) + h
Call subguoceng(x, y, xio)

For k = 1 To N
A(k, j) = y(k)
Next k
x(j) = x(j) - h
Next j

Call subguoceng(x, y, xio)

Call gaussj(A, y, N)

p = 1.0
For i1 = 1 To N
p = p - y(i1)
Next i1

For i2 = 1 To N
x(i2) = x(i2) - y(i2) * h / p
Next i2
h = h * w
Exit For
End If
Next i
Loop While yesno = False
Catch keyException As myException
pgkeyexception.Source = pgkeyexception.Source + "<错误讯息>：error:gyzcx.pnlfs" + Chr(10)
Throw pgkeyexception
End Try

End Sub

Public Sub gaussj(ByRef A(,) As Double, ByRef b() As Double, ByRef n As Integer)

Dim i, j, k, l, ll, irow, icol As Integer
Dim IPIV(n), INDXR(n), INDXC(n) As Integer

Dim big, dum, pivinv As Double

For j = 1 To n
IPIV(j) = 0
Next j
For i = 1 To n
big = 0.0#
For j = 1 To n
If IPIV(j) <> 1 Then
For k = 1 To n
If IPIV(k) = 0 Then
If System.Math.Abs(A(j, k)) >= big Then
big = System.Math.Abs(A(j, k))
irow = j
icol = k
End If
ElseIf IPIV(k) > 1 Then

pgkeyexception.Source = "<错误讯息>：error:gyzcx.gaussj" + Chr(10)
Throw pgkeyexception

Exit Sub
End If
Next k
End If
Next j
IPIV(icol) = IPIV(icol) + 1
If irow <> icol Then
For l = 1 To n
dum = A(irow, l)
A(irow, l) = A(icol, l)
A(icol, l) = dum
Next l
dum = b(irow)
b(irow) = b(icol)
b(icol) = dum
End If
INDXR(i) = irow
INDXC(i) = icol
If A(icol, icol) = 0.0# Then
pgkeyexception.Source = "<错误讯息>：error:gyzcx.gaussj" + Chr(10)
Throw pgkeyexception
Exit Sub
End If
pivinv = 1.0# / A(icol, icol)
A(icol, icol) = 1.0#
For l = 1 To n
A(icol, l) = A(icol, l) * pivinv
Next l
b(icol) = b(icol) * pivinv
For ll = 1 To n
If ll <> icol Then
dum = A(ll, icol)
A(ll, icol) = 0.0#
For l = 1 To n
A(ll, l) = A(ll, l) - A(icol, l) * dum
Next l
b(ll) = b(ll) - b(icol) * dum
End If
Next ll
Next i
For l = n To 1 Step -1
If INDXR(l) <> INDXC(l) Then
For k = 1 To n
dum = A(k, INDXR(l))
A(k, INDXR(l)) = A(k, INDXC(l))
A(k, INDXC(l)) = dum
Next k
End If
Next l

End Sub
Public Sub llrswd()
Dim llrswdsub As New gyzcx.pnlfsDelegate(AddressOf llrswdDelegate)
pnlfs(N, x, y, llrswdsub, xio)
end sub

Sub llrswdDelegate(ByRef x() As Double, ByRef y() As Double, ByRef xio() As Double)
'在这里写方程组代码
'y(1)=......
:
:
'y(n)=......

end sub

[Montaque]

看你用什么算法，比如高斯消元等等。形成倒三角。

[block]

同理可以使用下面的方程组代码
x0=x(1)
x1=x(2)
x2=x(3)
'其它已知量（a,b,c,d,e,n）在这里赋值
y(1)=n/b*x0+x1*exp(-x2)-a/b
y(2)=exp(-x2)*(x0+x1-c/b)
y(3)=x0*x1*exp(-x2)+x1^2*exp(-2*x2)-d/e*x1*exp(-x2)

顺便说明一下前面帖子中给出的那个代码中eps = Val(1 / (10 ^ rlmaina.jsjdszBar.Value))
你可以直接设置为eps =0.00001即可，Val(1 / (10 ^ rlmaina.jsjdszBar.Value))是原来程序中通过界面控制改变计算精度的，一般设置为0.00001即可。

[block]

将下面的方程组代码加入到前面的llrswdDelegate中即可
a=x(1)
b=x(2)
r=x(3)
'其它已知量在这里赋值
y(1)=∑yi-n*a-b∑exp(-r*xi)
y(2)=∑yi*exp(-r*xi)-a*∑exp(-r*xi)-b*∑exp(-2*r*xi)
y(3)=b*∑xi*yi*exp(-r*xi)-b*a∑exp(-r*xi)-b*b∑xi*EXP(-2*r*xi)

[Montaque]

不好意思，高斯消元只能处理线性的。
对于非线性，你可以采用其他的全局求优的算法，比如遗传算法，pso等等。

[ForMyLove]

整理后得：
n/b*x0+x1*exp(-x2)-a/b=0
exp(-x2)*(x0+x1-c/b)=0
x0*x1*exp(-x2)+x1^2*exp(-2*x2)-d/e*x1*exp(-x2)=0
现在的未知量是：
x0,x1,x2.
解出后，分数还可以再给。本人说话算话。

[ForMyLove]

方程组代码要怎么写呀？能是通用的吗？

[block]

方程组代码已经给你写出来了，就是2003-3-18 13:45:48 那帖。如果要通用就用前面的代码，使用托管，不同的方程组调用不同的托管过程。