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

_