牛顿迭代法求解线性方程组

[fyhdhr0022]

牛顿迭代法求解线性方程组的通用C算法，急需，谢谢大家了

[laolaoliu2002]

#include <stdio.h>
#include <math.h>
int Didai(g,x,eps,N)
double (*g)();
double *x,eps;
int N;
{
int i=0;
double x0=*x;
while (1)
{
*x=g(x0);
if (fabs(*x-x0)<eps) return 1;
if (++i>=N) return 0;
x0=*x;
}
}

double f1(double x)
{ return sqrt(10/(4+x)); }

double f2(double x)
{ return sqrt(10-x*x*x)/2; }

double f3(double x)
{ return x-(x*x*x+4*x*x-10)/(3*x*x+8*x); }

main()
{
double x=1.5,y=1.5,z=1.5;
if(Didai(f1,&x,1e-5,1000))
printf("root1=%f\n",x);
if(Didai(f2,&y,1e-5,1000))
printf("root2=%f\n",y);
if(Didai(f3,&z,1e-5,1000))
printf("root3=%f\n",z);
}

[laolaoliu2002]

#include <stdio.h>
#include <math.h>
int Didai(g,x,eps,N)
   double (*g)();
   double *x,eps;
   int N;
{
   int i=0;
   double x0=*x;
   while (1)
   {
      *x=g(x0);
      if (fabs(*x-x0)<eps) return 1;
      if (++i>=N) return 0;
      x0=*x;
   }
}
 
double f1(double x)
{ return sqrt(10/(4+x)); }
 
double f2(double x)
{ return sqrt(10-x*x*x)/2; }
 
double f3(double x)
{ return x-(x*x*x+4*x*x-10)/(3*x*x+8*x); }
 
main()
{
   double x=1.5,y=1.5,z=1.5;
   if(Didai(f1,&x,1e-5,1000))
      printf("root1=%f\n",x);
   if(Didai(f2,&y,1e-5,1000))
      printf("root2=%f\n",y);
   if(Didai(f3,&z,1e-5,1000))
      printf("root3=%f\n",z);
}

呵呵忘记贴成code形式的了,补一个

[laolaoliu2002]

 #include<iostream.h>   
  #include<math.h>   
    
  double   f1(double   x);       //声明原函数   
  double   f2(double   x);       //声明导函数   
    
  void   main()   
  {   
          double   a,x0;x1;   
    
          double   eps=1e-8;         //精度要求   
    
          cout<<"please   input   x0   :";   
          cin>>x0;   
    
          while(1)   
          {   
                  if(f2(x0)==0)             //计算x0处的导函数,若为0.则返回   
                  return;   
    
                  x1=x0-f1(x0)/f2(x0);   
                  if(abs(x1-x0)<eps)             //若满足精度要求就输出   
                  {   
                            cout<<"x1="<<x1<<","   
                                    <<"f1(x1)="<<f1(x1)<<","   
                                    <<endl;   
                  }   
                    break;                   //若满足精度要求就跳出       
                    elae                       //否则就继续循环   
                            x0=x1;   
              };   
  }   
    
  double   f1(double   x)             //定义原函数   
  {   
            //此处声明原函数体   
  }   
    
  double   f2(double   x)         //定义导函数   
  {   
            //此处声明导函数体   
  }                   

[fyhdhr0022]

如果是使用自己求导的函数的话我也有源程序，但是我现在是不明白如何使用C算法实现求导的运算，从而不通过自己的计算得到求导后的方程，谢谢大家

[fyhdhr0022]

新手上路，不明白如何在C中实现求导的运算。
如果是光通过在自己求导运算得出求导后的f(x)，那样的运算量也不是很少，对于一个高次的来说，也是很麻烦的。

[laolaoliu2002]

如果是求符号的话,很难!
比如只是多项式求导的话,就涉及多项式的表示方法(数据结构里有)
如果是其它函数的话,就更复杂了,复杂之处在于有什么样的数据结构来表示一个函数

[laolaoliu2002]

#define MAX   100
#define NULL 0
#define eps 0.00000001

double aa[MAX],bb[MAX];
int qj_num;

#include "math.h"
#include "conio.h"

typedef struct node
  {double coef;
   double exp;
   struct node *next;
   }NODE;


NODE *create()
{
  NODE *head,*r,*s;
  int coef,exp;
  clrscr();
  head=(NODE *)malloc(sizeof(NODE));
  r=head;
  printf("\nplease input coef(m) and exp(n):\n");
  scanf("%d%d",&coef,&exp);
  while(coef!=0)
  {
    s=(NODE*)malloc(sizeof(NODE));
    s->coef=coef;
    s->exp=exp;
	    r->next=s;
    r=s;
    printf("please input coef(m) and exp(n):\n");
    scanf("%d%d",&coef,&exp);
   }
  r->next=NULL;
  head=head->next;
  return(head);
}




NODE *create_qd(NODE *head)
{
  NODE *qdhead,*r,*s,*p;
  qdhead=(NODE *)malloc(sizeof(NODE));
  r=qdhead;
  p=head;
  for(;p!=NULL;p=p->next)
   {
     if(p->exp!=0)
      { s=(NODE *)malloc(sizeof(NODE));
	s->coef=p->coef*p->exp;
	s->exp=p->exp-1;
	r->next=s;
	r=s;
      }
     else
      break;
    }
   r->next=NULL;
   qdhead=qdhead->next;
  return(qdhead);
}


NODE *create_qdqd(NODE *qdhead)
{
  NODE *qdqdhead,*r,*s,*p;
  qdqdhead=(NODE *)malloc(sizeof(NODE));
  r=qdqdhead;
  p=qdhead;
  for(;p!=NULL;p=p->next)
   {
     if(p->exp!=0)
      { s=(NODE *)malloc(sizeof(NODE));
	s->coef=p->coef*p->exp;
	s->exp=p->exp-1;
	r->next=s;
	r=s;
      }
     else
      break;
    }
   r->next=NULL;
   qdqdhead=qdqdhead->next;
  return(qdqdhead);
}


double qz(NODE *p,double z)
{int i,j;
 double sum,n=1;
 for(j=0;j<p->exp;j++)
       n=n*z;
       sum=n*p->coef;
 return sum;
 }





double fx(x,head)
double x;
NODE *head;
{
  double fa=0.0;
  NODE *p;
  for(p=head;p!=NULL;p=p->next)
   fa=fa+qz(p,x);
	       return(fa);
}


void qujian(double a,double b,NODE *head)
{ void output_qujian();
  double  fa,fb;
  int flag=0,j=1;
   double i,w=1.0;
   double temp;


  while(flag==0&&j<=2)
   {   i=a;
       fa=fx(i,head);
       for(;i<=b;i=i+w)
	{
	      fb=fx(i,head);
	      if(fa*fb<=0)
		    {
		       if(fa==0&&fb==0)
			 {
			   fa=fb;
			   continue;
			 }
		       else
			 {
			   bb[qj_num]=i;
			   aa[qj_num++]=i-w;
			   flag=1;
			 }
		     }
	  fa=fb;
	}
   if(flag==1)
     printf("\nthis program have square root!");
  else
      {w=w/2;
       ++j;
       }
}
if(qj_num==0)
 printf("this ploy no square root!");
  else
   output_qujian();
}

void output_qujian()
{ int i;
 printf("\nqj_num is\n");
 for(i=0;i<qj_num;i++)
  {printf("[%lf %lf]\t",aa[i],bb[i]);
  }
}




double diedai_dian(double aa,double bb,NODE *head,NODE *qdqdhead)
{
   double a1=aa,b1=bb;
   double x;
   double half1,half2;
   double fa,fb;
   NODE *p=head;
   x=(a1+b1)/2;
   half1=fx(x,p);
   if(half1>0)
	 b1=x;
   if(half1<0)
	 a1=x;

   x=(a1+b1)/2;
   half2=fx(x,p);
   if(half2>0)
       b1=x;
   if(half2<0)
       a1=x;

 fa=fx(a1,head);
 fb=fx(a1,qdqdhead);
 if(fa*fb>0)
  return a1;
  else
   return b1;

 }



double newton(double xi,NODE *head,NODE *qdhead)
{
  double fxi,fxi1;
  double x;

  for(;;)
    {
      x=xi;
      fxi=fx(xi,head);
      fxi1=fx(xi,qdhead);
      xi=x-fxi/fxi1;
      if(fabs(xi-x)<=eps)
	 break;
     }
   return(x);
}

main()
{
 NODE *head=NULL,*qdhead=NULL,*qdqdhead=NULL;
  int i,a,b;
  double xi,x;
  clrscr();

  printf("please input a   b!");
  scanf("%d%d",&a,&b);
  head=create();
  qdhead=create_qd(head);
  qdqdhead=create_qdqd(qdhead);

  qujian(a,b,head);


  printf("\n");
  for(i=0;i<qj_num;i++)
  {  printf("[%lf %lf] \t",aa[i],bb[i]);
     xi=diedai_dian(aa[i],bb[i],head,qdqdhead);
     printf("%lf\t",xi);

     x=newton(xi,head,qdhead);
     printf("%lf\t",x);
     printf("\n");
  }

[laolaoliu2002]

呵呵，费好大劲