newton raphson算法

dyjun888666

哪位能提供一下newton raphson算法的程序阿，求解非线性方程组的算法，我在网上找不到。

风花雪月

一个matlab实例

1. % This m-file takes care of synthetic division.
   
2. % By giving one polynomial and one root this function returns
   
3. % the polynomial formed with the other roots of the given polynomial excluding the given root.
   
4. % Keerthi Venkateswara Rao-M.Tech-NITC-INDIA
   
5. function coeff_second=syn_division(coeff_function,fun_root_new)
   
6. order_fun=size((coeff_function),2);
   
7. coeff_second=0;
   
8. for index=1:size((coeff_function),2)-1
   
9.     if index==1
   
10.         coeff_second(index)=coeff_function(index);
    
11.     else
    
12.         coeff_second(index)=coeff_function(index)+fun_root_new*coeff_second(index-1);
    
13.     end
    
14. end

复制代码

1. % This m-file calculates the derivative of the function, the limitation of
   
2. % this function is, it can calculate only the derivatives of power(x,n)....
   
3. % Keerthi Venkateswara Rao-M.tech-NITC-INDIA.
   
4. function coeff_derivative=derivate(coeff_function)
   
5. der_order=size((coeff_function),2)-1;
   
6. coeff_derivative=0;
   
7. for index=1:size((coeff_function),2)-1
   
8.     coeff_derivative(index)=der_order*coeff_function(index);
   
9.     der_order=der_order-1;
   
10. end

复制代码

1. % this m-file calculates the real roots of the given polynomial using
   
2. % newton raphson technique.this m-file calls the functions in the two m-files named as syn_division and
   
3. % derivate.
   
4. % coeff_function is the 1xn matrix conatining the coeff of the polynomial.
   
5. % Keerthi Venkateswara Rao-M.Tech-NITC-INDIA
   
6. function [final_roots,functionvalue] = newton(coeff_function,initial_guess,error_tolerance,max_iterations)
   
7. iterations=0;
   
8. max_no_roots=size(coeff_function,2);
   
9. final_roots=0;
   
10. functionvalue=0;
    
11. for no_roots=1:max_no_roots-1
    
12.     fun_root_new=initial_guess;
    
13.     flag=1;
    
14.     coeff_der_function=derivate(coeff_function);
    
15.     order_fun=size(coeff_function,2);
    
16.     order_der_fun=size(coeff_der_function,2);
    
17.     while flag==1
    
18.         fun_root_old=fun_root_new;
    
19.         fx=0;
    
20.         dfx=0;
    
21.         nonzero=1;
    
22.         while nonzero==1
    
23.             powers=order_fun-1;
    
24.             for index=1:order_fun
    
25.                 fx=fx+coeff_function(index)*fun_root_old^powers;
    
26.                 powers=powers-1;
    
27.             end
    
28.             powers=order_der_fun-1;
    
29.             for index=1:order_der_fun
    
30.                 dfx=dfx+coeff_der_function(index)*fun_root_old^powers;
    
31.                 powers=powers-1;
    
32.             end
    
33.             if dfx==0
    
34.                 fun_root_old=fun_root_old+1;
    
35.             else
    
36.                 nonzero=0;               
    
37.             end               
    
38.         end
    
39.         iterations = iterations + 1;
    
40.         fun_root_new = fun_root_old - fx/dfx;
    
41.         if iterations >= max_iterations
    
42.             flag=0;
    
43.         elseif  abs(fun_root_new-fun_root_old)<=error_tolerance
    
44.             flag=0;
    
45.             final_roots(no_roots)=fun_root_new;
    
46.             functionvalue(no_roots)=fx;
    
47.         end
    
48.     end
    
49.     coeff_function=syn_division(coeff_function,fun_root_new);
    
50. end

复制代码

风花雪月

c语言实例

1. #include <math.h>
   
2. #include <iomanip.h>
   
3. #include <iostream.h>
   
4. #include <process.h>
   
5. 
   
6. void usrfun(double x[16],double alpha[16][16],double beta[16])
   
7. {
   
8.     int np;
   
9.     np = 15;
   
10.     alpha[1][1] = -2.0 * x[1];
    
11.     alpha[1][2] = -2.0 * x[2];
    
12.     alpha[1][3] = -2.0 * x[3];
    
13.     alpha[1][4] = 1.0;
    
14.     alpha[2][1] = 2.0 * x[1];
    
15.     alpha[2][2] = 2.0 * x[2];
    
16.     alpha[2][3] = 2.0 * x[3];
    
17.     alpha[2][4] = 2.0 * x[4];
    
18.     alpha[3][1] = 1.0;
    
19.     alpha[3][2] = -1.0;
    
20.     alpha[3][3] = 0.0;
    
21.     alpha[3][4] = 0.0;
    
22.     alpha[4][1] = 0.0;
    
23.     alpha[4][2] = 1.0;
    
24.     alpha[4][3] = -1.0;
    
25.     alpha[4][4] = 0.0;
    
26.     beta[1]= x[1]*x[1]+x[2]*x[2]+x[3]*x[3]-x[4];
    
27.     beta[2]=-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[4]*x[4]+1.0;
    
28.     beta[3]=-x[1]+x[2];
    
29.     beta[4]=-x[2]+x[3];
    
30. }
    
31. 
    
32. void main()
    
33. {
    
34.     //program d10r13
    
35.     //driver for routine mnewt
    
36.     int i,j,kk,k,ntrial,np,n;
    
37.     double tolx,tolf,xx;
    
38.     ntrial = 5;
    
39.     tolx = 0.000001;
    
40.     n = 4;
    
41.     tolf = 0.000001;
    
42.     np = 15;
    
43.     double x[16],alpha[16][16],beta[16];
    
44.     for (kk = -1; kk<= 1; kk=kk+2)
    
45.     {
    
46.         for (k =1; k<=3; k++)
    
47.         {
    
48.             xx = 0.2 * k * kk;
    
49.             cout<< "Starting vector number "<<k<<endl;
    
50.             for (i = 1; i<=4; i++)
    
51.             {
    
52.                 x[i] = xx + 0.2 * i;
    
53.                 cout<<"x("<<i<<")= "<<x[i]<<endl;
    
54.             }
    
55.             for (j = 1;j <=ntrial; j++)
    
56.             {
    
57.                 mnewt(1,x,n,tolx,tolf);
    
58.                 usrfun(x,alpha,beta);
    
59.                 cout<<"  i        x(i)            f"<<endl;
    
60.                 for (i = 1; i<=n; i++)
    
61.                 {
    
62.                     cout<< setw(3)<< i;
    
63.                     cout<< setw(14)<< x[i];
    
64.                     cout<< setw(16)<< -beta[i]<<endl;
    
65.                 }
    
66.             }
    
67.         }
    
68.     }
    
69. }

复制代码

1. void mnewt(int ntrial,double x[16],int n,double tolx,double tolf)
   
2. {
   
3.     double alpha[16][16],beta[16];
   
4.     int k,i,indx[16];
   
5.     double errf,d,errx;
   
6.     for (k = 1; k<=ntrial; k++)
   
7.     {
   
8.         usrfun(x,alpha,beta);
   
9.         errf = 0.0;
   
10.         for (i = 1; i<=n; i++)
    
11.             errf = errf + fabs(beta[i]);
    
12.         if (errf <= tolf)  _c_exit();
    
13.         ludcmp(alpha,n,indx,d);
    
14.         lubksb(alpha,n,indx,beta);
    
15.         errx = 0.0;
    
16.         for (i = 1; i<=n; i++)
    
17.         {
    
18.             errx = errx + fabs(beta[i]);
    
19.             x[i] = x[i] + beta[i];
    
20.         }
    
21.         if (errx <= tolx) _c_exit();
    
22.     }
    
23. }

复制代码

ChaChing

好资料!
2F的原始来源连接
http://www.mathworks.com/matlabcentral/fileexchange/4313

zjj04640451

谢谢分享，已经有些概念了

ME!

这也太复杂了吧，我的求解牛顿迭代的时候出现了复数，是初始值不对还是精度不对啊

redplum

的确复杂了，有些还是C++