506 哪位大哥哥能帮我把这个BASIC版的转成C语言版啊,救救小妹啊!!!! 用最速下降法解非线性方程组
一. 功能
求非线性实函数方程组
fj(x1,x2,…,xn)=0,j=1,2,…,n (1)
的若干组实数解,因fj(x1,x2,…,xn)=0表示曲线或曲面,一组解
就是它们的一个交点,当然交点往往不止一个.
二. 算法简介
该方法的基本思想,是平方和函数
F(x1,x2,…,xn)=∑fj^2(x1,x2,…,xn) (2)
的一个实零点.若F(x1,x2,…xn)<E,则认为该x1,x2,…,xn是
方程组(1)的满足精度要求的一组解.它可看作F(x1,x2,…,xn)
的极小值点.E可取作E=1E-5或1E-6等
具体计算步骤
1. 给出一组初值 X^(0)=[X1^(0),X2^(0),…,Xn^(0)]
其中Xj^(0) (j=1,2,…,n)不能全为0
2. 从点X^(0)=[X1^(0),X2^(0),…,Xn^(0)]出发,沿其梯度
方向的相反方向去寻找函数F(x1,x2,…,xn)的极小值,设
X^(0)处的梯度为
G=(ЭF/Эx1,ЭF/Эx2,…,ЭF/xnЭ)^T (3)
其中ЭF/Эxi={F[x1^(0),…,xi-1^(0)xi^(0) +hi,
xi+1^(0),…,xn^(0)]-F[x1^(0),x2^(0),…,xn^(0)]}/
hi,I=1,2,…,n (4)
在(4)式中可取
hi=Cxi^(0),xi^(0)≠0
hi=C,xi^(0)=0
即C可选为10^-5----10^-1,不能太小.
3. 计算x^(1)=[x1^(1),x2^(1),…,xn^(1)]
xi^(1)=λЭF/Эxi,(I=1,2,…,n)
其中λ1计算比较复杂,它确定着沿X^(0)处梯度负方向以
多大步长向前跨进,可用下式近似代替
λ1=F[x1^(0),x2^(0),…,xn^(0)]/∑(ЭF/Эxi),
I=1,2,…,n
4. 计算F[x1^(1),x2^(1),…,xn^(1)],若<E,则X^(1)为所
求,否则再从x^(1)=[x1^(1),…,xn^(1)]出发以该点梯度
的反方向向前寻找,即以点X^(1)代替点X^(0),重复2-4的手续.已迭代k次,F[x1^(k),…,xn^(k)]<E,则x^(k)=
[x1^(k),…,xn^(k)]就是与所给初值X^(0)相应的解.
注:当初值不一样时,可能得到不同的解.要求方程组的若干组解,就得给出不同的初值进行计算,该方法对初值都是收敛的,但收敛速度不同,且都是线性的,即比较慢.
三. 程序使用说明
1. 输入参数
N 方程组中方程的个数(也是变元的个数)
E 精度控制常数,可取E=1E-5---1E-10
H 精度控制常数,求数值偏导数时个变量的增量,可取h=
1E-5---1E-7,不能太小,程序有时需用它作为除数,太小会出
现溢出的现象.
Y 程序中于330---360语句中分别放着两个方程组的求解初值,当Y=1时求第一个方程的解,当Y=2时求第二个方程的解.
2. 输出参数
初值,方程组及其解,均用汉字说明.
3. 在求解第二个方程组时,运行前需消区该方程初值(存于360
语句中)前面存放数据的330,340语句行.
四.结果分析
例1,
输入数据:0.1,0.5,–0.4
输出结果:2.1157133485E-4,1.5492098024,-2.78187035E-3
示例结果:4.266791E-5,1.549253,4.96848E-4
例2,
输入数据:1.2 ,6.1,-5
输出结果:1.0001580128,5.0003882579,-4.0003406799
示例结果: 1.000163,5.000381,-4.000338
例3,
输入数据:0.1,0.5
输出结果:2.3271791720E-1,5.6610944001E-2
示例结果: 0.2326433,5.651056E-2
分析:结果与示例有偏差,其数量级为E-4,与精度的选择有关,其结论应为正确。
Basic程序:
10'************************
20'506
30'************************
40 L=1
50 INPUT "n=";N:INPUT "E,h=";E,H
60 INPUT "Y=";Y
70 DIM X(N),H(N),D(N),Y(N)
80 PRINT TAB(3);"FIRST NUM:"
90 FOR I=1 TO N
100 READ X(I):PRINT X(I);
110 NEXT I
120 IF L=1 THEN 130 ELSE 280
130 IF Y<>1 THEN 220
140 DEF FNA(X)=X(1)-5*X(2)^2+7*X(3)^2+12
150 DEF FNC(X)=3*X(1)*X(2)+X(1)*X(3)-11*X(1)^2
160 DEF FNB(X)=2*X(2)*X(3)+40*X(1)
170 PRINT TAB(3);"EQUATIONS:"
180 PRINT "fnb(x)=x(1)-5*x(2)^2+7*x(3)^2+12"
190 PRINT "fnb(x)=3*x(1)*x(2)+x(1)*x(3)-11*x(1)^2"
200 PRINT "fnc(x)=2*x(2)*x(3)+40*x(1)"
210 GOTO 280
220 IF L=1 THEN 230 ELSE 280
230 DEF FNAA(X)=4*X(1)-X(2)+EXP(X(1))/10-1
240 DEF FNBB(X)=-X(1)+4*X(2)+X(1)^2/8
250 PRINT TAB(3);"EQUATIONS:"
260 PRINT "fnaa(x)=4*x(1)-x(2)+exp(x(1))/10-1"
270 PRINT "fnbb(x)=-x(1)+4*x(2)*x(1)^2/8"
280 GOSUB 500
290 PRINT TAB(3);"JIE:"
300 FOR I=1 TO N
310 PRINT X(I);
320 NEXT I
330 DATA 0.1,0.5,-0.4
340 DATA 1.2,6.1,-5.0
350 '*****
360 DATA 0.2,0.05
370 IF Y=2 THEN 400
380 L=L+1
390 IF L=2 THEN 80
400 END
410'***
420 IF Y<>1 THEN 450
430 Y(1)=FNA(X):Y(2)=FNB(X):Y(3)=FNC(X)
440 RETURN
450 Y(1)=FNAA(X):Y(2)=FNBB(X)
460 RETURN
500 'ZI CHENG XUE'
510 FOR I=1 TO N
520 IF X(I)<>0 THEN 550
530 H(I)=H
540 GOTO 560
550 H(I)=H*X(I)
560 NEXT I
570 GOSUB 410
580 F=0
590 FOR J=1 TO N
600 F=F+Y(J)^2
610 NEXT J
620 IF (F-E)<=0 THEN 800
630 S=0
640 FOR I=1 TO N
650 X(I)=X(I)+H(I)
660 GOSUB 410
670 F0=0
680 FOR J=1 TO N
690 F0=F0+Y(J)^2
700 NEXT J
710 D(I)=(F0-F)/H(I)
720 S=S+D(I)^2
730 X(I)=X(I)-H(I)
740 NEXT I
750 T=F/S
760 FOR I=1 TO N
770 X(I)=X(I)-T*D(I)
780 NEXT I
790 GOTO 510
800 RETURN
_