求助！！ 数值分析程序（c语言）

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
<>1。牛顿迭代法</P>
  d. W- X, Y8 }: ]0 k' y; X<>1。全主消元法</P>
<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
( H& z* ~" D! P$ d9 w<>拜托了。。。。</P>

loveyuan

<>1。牛顿迭代法</P><>
  #include "stdio.h"
- E; t+ `3 Z) `; ^: y+ }8 a; O  #include "math.h"
4 R* c6 ~& \2 p) z/ i/ j9 l+ C  int dnewt(x,eps,js)
  int js;
8 d9 v( [, s5 [# c3 f0 t  double *x,eps;
  { extern void dnewtf();
3 U6 N8 h! G1 T4 [  ]/ u% J    int k,l;
" Z! _$ f( |/ |6 v/ q0 J+ m    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
! p* x# A  r8 @; `) y/ K1 E    dnewtf(x0,y);
    d=eps+1.0;
    while ((d&gt;=eps)&amp;&amp;(l!=0))
      { if (fabs(y[1])+1.0==1.0)
+ U' u2 q, k1 _          { printf("err\n"); return(-1);}
        x1=x0-y[0]/y[1];
/ L- B: r5 b& |9 h( T        dnewtf(x1,y);
9 {. x5 \( M6 @) N% E; a        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
      }
    *x=x1;
    k=js-l;
    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
3 h+ z# D; |; s. d  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
0 L* L) Z8 R2 L& ^  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
0 q9 J* u7 c! Y2 H0 ^2 n' l    js=malloc(n*sizeof(int));
    for (k=0;k&lt;=n-2;k++)
+ W" W. n5 o, Z$ d      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
        for (j=k;j&lt;=n-1;j++)
          { u=i*n+j;
            p=ar*ar+ai*ai;
            if (p&gt;d) {d=p;js[k]=j;is=i;}
          }
        if (d+1.0==1.0)
5 m! ?. o  x( D          { free(js); printf("err**fail\n");
* |5 L" @* \8 c7 g/ j            return(0);
8 [# Q+ Z1 ?; M/ q% A! q          }
& @7 @9 h) M. z8 D7 V2 K        if (is!=k)
1 K4 q3 D3 L, O& H3 ~# D          { for (j=k;j&lt;=n-1;j++)
+ q9 }. h* j( A+ g# d7 J) _              { u=k*n+j; v=is*n+j;
                p=ar; ar=ar[v]; ar[v]=p;
                p=ai; ai=ai[v]; ai[v]=p;
              }
            p=br[k]; br[k]=br[is]; br[is]=p;
. k$ s( @: B0 j5 n            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
5 i% W& T' `3 y/ T          }
3 m: F$ j, G- a. j+ g        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
3 L3 n9 e4 w2 j6 t! u+ y1 e! d1 g              p=ar; ar=ar[v]; ar[v]=p;
              p=ai; ai=ai[v]; ai[v]=p;
            }
        v=k*n+k;
1 N- ?. c  ~0 J        for (j=k+1;j&lt;=n-1;j++)
: a1 g/ ?$ `: t0 T$ V          { u=k*n+j;
& x. C  ?1 N% t, ~+ k, ^2 D: A            p=ar*ar[v]; q=-ai*ai[v];
; A: X- \& q# b7 H            s=(ar[v]-ai[v])*(ar+ai);
) R9 m; W+ h' j' S# v            ar=(p-q)/d; ai=(s-p-q)/d;
+ |% Z3 k+ q3 j: F          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
& b& @: m) f: c        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i&lt;=n-1;i++)
6 a# B; Y, [4 A1 W' e: q          { u=i*n+k;
# t# b5 q9 t( K1 ]: q" D  V            for (j=k+1;j&lt;=n-1;j++)
              { v=k*n+j; l=i*n+j;
( G7 E' l1 ]2 T" l                p=ar*ar[v]; q=ai*ai[v];
                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
( }- H* j8 X* K: |5 O* K                ai[l]=ai[l]-s+p+q;
  E& v( n3 h8 p7 `* S- s              }
" k4 n) e- O& ~, R            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
- x3 e0 k5 e% A+ y7 L            br=br-p+q; bi=bi-s+p+q;
          }
9 D% c9 [# n6 l      }
# F5 q, P+ f. G. a8 e    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
' H+ v- ]) }: ?6 ?' l    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
        return(0);
      }
    p=ar*br[n-1]; q=-ai*bi[n-1];
0 j$ G7 E0 ^* B/ Z5 S% ~! S    s=(ar-ai)*(br[n-1]+bi[n-1]);
    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
    for (i=n-2;i&gt;=0;i--)
; s  i3 |: Q! r9 G% l    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
        p=ar*br[j]; q=ai*bi[j];
        s=(ar+ai)*(br[j]+bi[j]);
        br=br-p+q;
        bi=bi-s+p+q;
      }
" r* H. Y+ _6 i& \3 K# Z/ V. v    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
# y& h; _$ o! l$ N9 v1 |      if (js[k]!=k)
        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
: ^+ e' w( J4 V( q  ^        }
    free(js);
4 F  i0 C* {% Y( A( A    return(1);
  }</P><>平方根法</P><>#include "math.h"
: S$ N2 _" X# j  #include "stdio.h"
  int achol(a,n,m,d)
2 t2 s$ a4 ^" Z' p: [0 M4 |. B; }  int n,m;
* a* |3 J3 O2 M2 g% R0 y  double a[],d[];
  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
      { printf("fail\n"); return(-2);}
& z. J  J: q) A& \1 w    a[0]=sqrt(a[0]);
! W8 b7 t( f' v2 L7 E    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
4 ^( y4 l: [5 u    for (i=1; i&lt;=n-1; i++)
      { u=i*n+i;
( l. w0 j8 _1 \) \  d6 \        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
4 J, s" L* @0 Y- a& q2 _          }
        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
        a=sqrt(a);
        if (i!=(n-1))
6 R. g) D5 O1 B, S" z8 N, }, B8 c# h          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
! m( o) p' i" h                for (k=1; k&lt;=i; k++)
                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
                a[v]=a[v]/a;
              }
/ b8 @' y+ t" q" a: Z# O) g          }
      }
* A0 ?# w/ v/ L9 t0 B    for (j=0; j&lt;=m-1; j++)
) j) f) o1 S; c      { d[j]=d[j]/a[0];
5 ]" x1 r; n/ _* x" D        for (i=1; i&lt;=n-1; i++)
2 S2 w. ]) q' T7 {& Q          { u=i*n+i; v=i*m+j;
            for (k=1; k&lt;=i; k++)
              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
( n- I- G: {4 |' q! @& |            d[v]=d[v]/a;
          }
6 K+ N( B" }: W5 Y( V6 d0 [      }
/ @* E( B0 g9 F! T* h9 H( ?/ P) B    for (j=0; j&lt;=m-1; j++)
      { u=(n-1)*m+j;
        d=d/a[n*n-1];
        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
% M% S9 U6 r. g- m" Z            for (i=k; i&lt;=n-1; i++)
              { v=(k-1)*n+i;
. S9 v9 J+ ]/ ~# i# `                d=d-a[v]*d[i*m+j];
+ Y" T+ I) w+ x* L4 w3 E1 a% f              }
            v=(k-1)*n+k-1;
            d=d/a[v];
          }
      }
# G; W+ b- p1 ^) a) @6 ]    return(2);
; t( q3 O1 @. L+ I. |6 I  Z  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
+ K9 y' v# p. @! D1 T  int n;
  ?6 e; z! o, Z; T: H# V% o  double x0,h,t,y[];
  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
    z=0.0;
3 v+ l. o$ z1 e; U    if (n&lt;1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return(z);
      }
    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
6 K- n7 m* o+ O      }
    else i=0;
6 }8 X* m! e1 r" G, x5 n    k=i-4;
& `/ u9 E3 O* h    if (k&lt;0) k=0;
    m=i+3;
    if (m&gt;n-1) m=n-1;
    for (i=k;i&lt;=m;i++)
      { s=1.0; xi=x0+i*h;
9 F$ G; E8 w9 V  B( m+ O# P        for (j=k; j&lt;=m; j++)
# V: ^0 d1 F1 w3 K# }# ?5 L/ [2 w4 N          if (j!=i)
            { xj=x0+j*h;
3 P8 _/ a+ @- l9 H1 \; ~5 P6 \              s=s*(t-xj)/(xi-xj);
' a. ?3 O1 U: b8 b$ g% Q+ {            }
* q  ?: O& S4 h        z=z+s*y;
      }
6 F+ B0 m  m4 M' U    return(z);
; v( x6 S; z- C( u+ i/ f  }
4 s* r% E4 O1 F. D3 ]1 m向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
3 G$ M: g9 F$ I) V) g  int n,m;
  double x[],y[],a[],dt[];
. N. [' c  y" O3 S  { int i,j,k;
& D: ?1 v) O9 @9 l* i" e# r    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
    for (i=0; i&lt;=m-1; i++) a=0.0;
    if (m&gt;n) m=n;
    if (m&gt;20) m=20;
    z=0.0;
    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
, [+ _# T3 u( G    for (i=0; i&lt;=n-1; i++)
      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
    a[0]=c*b[0];
    if (m&gt;1)
      { t[1]=1.0; t[0]=-p;
4 m2 l* e; U  R9 E6 R" u        d2=0.0; c=0.0; g=0.0;
" R3 R6 V9 E( D8 r( ]( d: M        for (i=0; i&lt;=n-1; i++)
          { q=x-z-p; d2=d2+q*q;
            c=c+y*q;
            g=g+(x-z)*q*q;
& V* u( e& k6 ^2 x7 a8 {          }
        c=c/d2; p=g/d2; q=d2/d1;
# J' t$ X6 k0 @( s$ _% |- @        d1=d2;
9 ]1 i- X, ]  s) h7 @- p/ B) m* v        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
. `$ V/ D1 K. P4 g# N  f    for (j=2; j&lt;=m-1; j++)
3 \1 t( R" g! d/ b+ i/ i  |& i. ~8 U      { s[j]=t[j-1];
5 C# v* p0 Y3 H        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
          for (k=j-2; k&gt;=1; k--)
  H6 P. k' o  m+ e1 o5 R* F            s[k]=-p*t[k]+t[k-1]-q*b[k];
        s[0]=-p*t[0]-q*b[0];
6 u( M9 N' b* ], P        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
4 ]1 D% n  _, X$ ]# r          { q=s[j];
            for (k=j-1; k&gt;=0; k--)
              q=q*(x-z)+s[k];
% _/ X. X* _1 @& f            d2=d2+q*q; c=c+y*q;
) i& ?1 f% D9 j* w            g=g+(x-z)*q*q;
( X/ l4 G# S, Q4 {( e          }
4 {4 `" A; \4 H- a' C        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
& G( \! W8 B# L6 P: j( _% V$ v        a[j]=c*s[j]; t[j]=s[j];
0 m9 h9 T+ C( g8 v% J" W        for (k=j-1; k&gt;=0; k--)
* I0 s$ d* d& x1 a! G4 W4 [4 j9 F          { a[k]=c*s[k]+a[k];
( h& e, x6 \, E7 p, E( F/ r! ~            b[k]=t[k]; t[k]=s[k];
! P! u# x' W0 m& v. H: E5 w8 b          }
+ B4 x& P( E7 u9 i1 B9 \+ O      }
2 U; `) h. k* w1 @7 p2 g( y) W    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
; L0 b; w6 v1 _7 F    for (i=0; i&lt;=n-1; i++)
6 X1 z/ X+ ^7 q3 Y      { q=a[m-1];
; Z% ^4 ^: `( g  s- n        for (k=m-2; k&gt;=0; k--)
- E% ]8 X2 N/ a8 {% D          q=a[k]+q*(x-z);
8 B9 z& A' M* J& b8 o  L" K. ~        p=q-y;
2 x& n" O. F- B  A9 |: i2 g        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
' ], U' y1 A% d! e4 c5 E1 [        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
      }
    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
  double fromb(a,b,eps)
6 H5 z" r+ v7 c' i1 K  double a,b,eps;
! |7 D6 h7 b# |0 Y3 u' `  { extern double frombf();
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
5 k) Y7 u6 V& U" `  E7 B    h=b-a;
0 {0 `  C( y! E6 D3 C    y[0]=h*(frombf(a)+frombf(b))/2.0;
- h% h( k7 v$ k    m=1; n=1; ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
      { p=0.0;
. N1 E2 G% v* ~        for (i=0;i&lt;=n-1;i++)
! P  Z) l3 v; }0 J( R( L& a          { x=a+(i+0.5)*h;
7 o- y; q4 v( @+ g            p=p+frombf(x);
, X9 m) w1 [! o8 g4 B' y          }
+ }! \: i5 j8 S  M) ?- w( V        p=(y[0]+h*p)/2.0;
# c; Z3 k4 W$ O$ ]        s=1.0;
7 y) {2 L9 M' f: R% B2 k& w$ p        for (k=1;k&lt;=m;k++)
! O8 Q- b, s8 ^2 e& m9 f1 s+ f# s          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
- ^# a# L9 _5 L  N$ B            y[k-1]=p; p=q;
: o" T" E- @! T) [' m          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
. Y1 h2 d  ~% C" n3 G/ ?- \      }
5 X: z( O$ q6 c6 C* ^$ s. I( |% z    return(q);
  }</P><>呵呵 希望对你有用！！</P>

西西

<>谢谢你哦，呵呵</P>

陈了

哇哇哇哇哇哇。。。好复杂呀！看不懂