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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
& p! S  B( [, A) T4 Q<>1。牛顿迭代法</P>
6 \1 z6 e& w* ]: ~- t<>1。全主消元法</P>
7 r) e, K( S8 m3 [5 N$ g' D<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
6 Z4 s* }9 {  y0 S9 Z<>6。龙贝格积分法</P>
<>拜托了。。。。</P>

陈了

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

西西

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

loveyuan

<>1。牛顿迭代法</P><>
0 c. J" h) ]. O' m5 f+ [5 o* D, Y% e  #include "stdio.h"
  z: _4 r  h2 C; C' n  #include "math.h"
8 R7 E( T4 ~& |1 B$ h1 Z' W2 T  int dnewt(x,eps,js)
  int js;
( T' N3 i0 _* ~( u4 ~+ r8 B  double *x,eps;
  { extern void dnewtf();
    int k,l;
8 e: F% b) \" @2 R4 V4 t! f    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
  }( _# k2 R0 j- T( J+ L9 @    dnewtf(x0,y);
2 W- W. r5 O/ o6 L' J    d=eps+1.0;
    while ((d&gt;=eps)&amp;&amp;(l!=0))
      { if (fabs(y[1])+1.0==1.0)
          { printf("err\n"); return(-1);}
6 E) u$ M0 e7 j" x        x1=x0-y[0]/y[1];
        dnewtf(x1,y);
        d=fabs(x1-x0); p=fabs(y[0]);
) n( p( o! N5 ^6 T        if (p&gt;d) d=p;
) d! l5 d# J; G        x0=x1; l=l-1;
' E3 J0 s. \/ _/ S  a' J" D      }
0 M: V- d2 g( i$ @/ N! G    *x=x1;
  G* F+ E+ O+ I3 j8 g    k=js-l;
    return(k);
/ n7 ?- l1 u4 E6 Z5 f  }</P><>全主消元法</P><>#include "stdlib.h"
2 x; E& O, X2 q. S( A" @; P  #include "stdio.h"
2 u1 ^3 b- i- W8 o" `  int acgas(ar,ai,n,br,bi)
& T; T: R- a. s  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
7 a- N6 Y& h( U' @& _: ]    for (k=0;k&lt;=n-2;k++)
) ~/ Z) ~7 n& p/ O* C      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
        for (j=k;j&lt;=n-1;j++)
! p. \- j* ~4 P1 k& D          { u=i*n+j;
0 ^! A% Q' T* r$ z            p=ar*ar+ai*ai;
            if (p&gt;d) {d=p;js[k]=j;is=i;}
          }
        if (d+1.0==1.0)
) b/ t7 V) V! x" l6 Q          { free(js); printf("err**fail\n");
" U! V5 Z# t3 z' v# T  v9 x. y            return(0);
, u. a# Y0 q* Y          }
/ E1 w/ ~& O+ T/ o6 R) u0 {        if (is!=k)
% [4 ~: ]) L% ^/ I2 A4 g5 w& h+ }) m          { for (j=k;j&lt;=n-1;j++)
, X  x& p3 J8 t; P$ W) a              { u=k*n+j; v=is*n+j;
" A; R( ?' ]5 }; N$ M4 c                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;
8 w4 E$ W1 S9 h2 Y: u            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
          }
        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
9 C$ W3 k2 Q: w2 T- f8 p- j              p=ar; ar=ar[v]; ar[v]=p;
- U" O" Z/ ]/ A              p=ai; ai=ai[v]; ai[v]=p;
% _) y8 _' ~$ F" K6 u* [! S' x            }
& Z1 N3 ^  H* p& x; v  [        v=k*n+k;
5 K: P8 I6 ~4 F1 j: v- k        for (j=k+1;j&lt;=n-1;j++)
9 W- ^* ^9 C0 y- M1 z  }8 X          { u=k*n+j;
2 J0 U3 g0 j9 p* Y9 Z5 v            p=ar*ar[v]; q=-ai*ai[v];
            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
          }
. A( p7 K! u& K( ^        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
0 E* ?3 m7 C. b* x# d' G% y        for (i=k+1;i&lt;=n-1;i++)
& P& o9 a& X( q9 U2 b          { u=i*n+k;
            for (j=k+1;j&lt;=n-1;j++)
+ C% X% _7 ?$ c              { v=k*n+j; l=i*n+j;
" q+ c0 K  U: S$ O1 Z4 l                p=ar*ar[v]; q=ai*ai[v];
! T& I- }" [; b6 R1 |$ K+ S% V                s=(ar+ai)*(ar[v]+ai[v]);
4 f/ @- s4 t, p                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
            br=br-p+q; bi=bi-s+p+q;
* k. [( H/ v& f, B0 X          }
      }
    u=(n-1)*n+n-1;
, E7 N  Z# i4 I) N6 @    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
# G7 s* M& c, A! v4 R1 q        return(0);
      }
* e4 t" W1 N' d" D4 z2 y9 T5 H    p=ar*br[n-1]; q=-ai*bi[n-1];
    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--)
    for (j=i+1;j&lt;=n-1;j++)
% w2 u7 K3 ^( p5 L      { u=i*n+j;
        p=ar*br[j]; q=ai*bi[j];
, S+ R5 q" p* M# ~3 T, A- n        s=(ar+ai)*(br[j]+bi[j]);
        br=br-p+q;
* N, V5 h* j! U        bi=bi-s+p+q;
! y1 e; a- K: O) i. A6 x8 P      }
/ O0 v; S+ |; \4 ~- p1 p1 x    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
5 w; J  J9 ]& G5 N3 [" A1 t      if (js[k]!=k)
+ R/ m! P# N: Q7 F        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
6 h! ^2 I  g, k8 d+ q5 m6 V        }
5 `2 }* w0 W: l    free(js);
    return(1);
  }</P><>平方根法</P><>#include "math.h"
, \- v/ n1 h) w7 x5 H( H$ g8 e  #include "stdio.h"
6 b" |, j) C* E1 O  int achol(a,n,m,d)
  int n,m;
' Y$ k: m. s5 W, p  double a[],d[];
) ~7 L# O4 \5 N2 C  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
' m! i; E+ u. o) X8 n* F      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
  n9 z5 Y# W4 ]4 X% {    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
- A/ C1 o# g0 f2 o/ l    for (i=1; i&lt;=n-1; i++)
' L9 l! E/ T- E! P      { u=i*n+i;
2 x. t; u8 j8 n: s5 O4 t+ e% w- L        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
, f8 Z. a; p* f$ u# N          }
; H/ S3 K% h5 d- _. Y' V8 T  A        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
        a=sqrt(a);
        if (i!=(n-1))
          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
                for (k=1; k&lt;=i; k++)
3 J' R% R; ?: S9 U) M2 ?                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
; t- e3 E: Q4 F! W. R: }                a[v]=a[v]/a;
              }
          }
      }
* r0 @6 n. a( R    for (j=0; j&lt;=m-1; j++)
/ g* y6 l6 _: ?. ^. @1 O8 B      { d[j]=d[j]/a[0];
* D' }: R" S4 c0 P        for (i=1; i&lt;=n-1; i++)
- p) l3 s2 Q3 [6 {( |          { u=i*n+i; v=i*m+j;
0 W& }: N6 U( y; v4 `            for (k=1; k&lt;=i; k++)
              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a;
, K, I7 O* J6 z. R          }
9 e9 h3 g0 p' B- I      }
+ F- L) c0 {6 ?- m" u4 g    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--)
; p$ W3 r. K" L% `$ z. o          { u=(k-1)*m+j;
            for (i=k; i&lt;=n-1; i++)
              { v=(k-1)*n+i;
- u! c" h# _1 D0 d* ~3 U) d& b6 O! m5 q" O                d=d-a[v]*d[i*m+j];
1 ]0 D/ c* c- v8 ^              }
            v=(k-1)*n+k-1;
            d=d/a[v];
6 A/ Z  V+ s3 A( [, |4 ]2 s          }
- F/ ]0 F( G$ c1 o. ~( {7 M9 A      }
3 X1 H# }8 J4 M6 F  s: o3 i( F7 Q    return(2);
  R% {' f5 s* A  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
' C0 t- K2 `3 W! X7 [0 l+ _* z  int n;
  double x0,h,t,y[];
7 e$ m' T$ n. ~0 B5 y* m  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
. N- m( d* L& ~; f1 {$ M    z=0.0;
7 N% _$ P/ l# {" C# d/ Z5 y7 [3 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)
3 ~1 K& [' c/ Q, f7 X      { p=(t-x0)/h; i=(int)p; q=(float)i;
: b4 L/ h; U  t# |        if (p&gt;q) i=i+1;
  z; y+ t/ R  M# E: V      }
/ x4 ?1 c+ B: P0 T    else i=0;
    k=i-4;
    if (k&lt;0) k=0;
    m=i+3;
# W$ T9 l' I$ w2 Q, ?    if (m&gt;n-1) m=n-1;
9 W; Y, c  ?4 y    for (i=k;i&lt;=m;i++)
/ R* p, H. \+ ]      { s=1.0; xi=x0+i*h;
) G7 M! A+ N! [; m" d3 }  E        for (j=k; j&lt;=m; j++)
4 c* O- V$ A- `% d, ^0 l          if (j!=i)
) r, _7 |' j, b+ F4 A0 Z, ^( g" j& B            { xj=x0+j*h;
              s=s*(t-xj)/(xi-xj);
8 N' |0 K, z- i; I% }            }
        z=z+s*y;
3 ]5 V, W% i% V# p7 [, J      }
    return(z);
  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
8 w) I$ R6 e0 n6 l. v6 z) r, j  int n,m;
* t, q/ Z/ ?* I/ z' I( D* [  double x[],y[],a[],dt[];
6 M5 W2 z9 P7 W7 h4 t  { int i,j,k;
4 ^3 p% q5 z% W* ], P    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;
% g. b- \3 l) G; w, y  T& v; d    if (m&gt;n) m=n;
    if (m&gt;20) m=20;
* y4 b! x& j, L: y( a( _' C  M    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;
    for (i=0; i&lt;=n-1; i++)
0 O* q. S8 @5 E+ Z. m4 Y. x# R3 x7 @* o      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
    a[0]=c*b[0];
  v$ f  f! O% E) P9 o& e7 s    if (m&gt;1)
      { t[1]=1.0; t[0]=-p;
/ g" n  d' o" s8 y        d2=0.0; c=0.0; g=0.0;
. `+ A4 [# h4 \$ |' F        for (i=0; i&lt;=n-1; i++)
          { q=x-z-p; d2=d2+q*q;
( o$ X  ^$ @7 q/ W8 i- T1 A2 t            c=c+y*q;
, G: Y' v7 ]1 m8 C! C            g=g+(x-z)*q*q;
          }
$ }+ Y5 ~0 {" ?        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
5 t! e" F3 ^: [4 J& d4 L        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
0 n$ ]0 V! G. ]' o" C9 P7 T+ K    for (j=2; j&lt;=m-1; j++)
  ~- J) y5 b$ V7 I      { s[j]=t[j-1];
0 j4 @1 M+ G- Y6 f0 R        s[j-1]=-p*t[j-1]+t[j-2];
- z8 G& Z' O0 k, f; r3 A3 O6 }5 a        if (j&gt;=3)
! \7 s$ W' r/ g          for (k=j-2; k&gt;=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
' v) z/ k4 Q+ U2 W- y        s[0]=-p*t[0]-q*b[0];
        d2=0.0; c=0.0; g=0.0;
3 [) [# b8 `' j3 @& {. U* G3 v        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
4 D# z# ~+ r- z' l  N, r  k            for (k=j-1; k&gt;=0; k--)
$ }; v3 S9 j4 c              q=q*(x-z)+s[k];
            d2=d2+q*q; c=c+y*q;
            g=g+(x-z)*q*q;
          }
        c=c/d2; p=g/d2; q=d2/d1;
1 C6 K7 i" x; t, E( K        d1=d2;
        a[j]=c*s[j]; t[j]=s[j];
; K: V" I% F1 V* U+ w" o        for (k=j-1; k&gt;=0; k--)
          { a[k]=c*s[k]+a[k];
8 t) u& g2 j% z* o; O2 \            b[k]=t[k]; t[k]=s[k];
          }
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
/ {( k2 c9 `# A/ ?    for (i=0; i&lt;=n-1; i++)
" ~1 ^: T* z. v      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
          q=a[k]+q*(x-z);
: k. r3 l/ [! ?0 R        p=q-y;
. c7 u' V4 w4 C7 g: l1 V        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
. w. x) T7 e7 Q! v$ o! O- N4 V        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
) W: m# i) X8 Q; y0 c# i      }
    return;
; F2 N$ |! x2 |7 ?% A  }</P><>龙贝格积分法</P><>#include "math.h"
+ z0 \' v+ h8 ?, ]/ y/ Q3 D6 J  double fromb(a,b,eps)
  double a,b,eps;
5 s! a$ n$ P) z8 b" x  { extern double frombf();
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
9 x! `: E$ ^+ ?3 U7 u* j" p    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
7 y) ?  ]. Q% k' r      { p=0.0;
% {* D" T9 H& j6 B  A$ i* g        for (i=0;i&lt;=n-1;i++)
          { x=a+(i+0.5)*h;
            p=p+frombf(x);
          }
* u& c/ W9 R3 f/ I/ c; {        p=(y[0]+h*p)/2.0;
* v; D# k8 I' j, M+ R3 c        s=1.0;
        for (k=1;k&lt;=m;k++)
          { s=4.0*s;
+ _6 t: B0 j' y3 A+ ^            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
2 u* v; l& F3 L! b. q          }
        ep=fabs(q-y[m-1]);
. `9 `% D3 ~7 G- o: R3 e3 E5 [        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
" C# T2 u! Z5 ], [      }
    return(q);
4 n& i2 \- H8 G$ F& o% @  }</P><>呵呵 希望对你有用！！</P>