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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
<>1。牛顿迭代法</P>
! m# L' J6 n' \<>1。全主消元法</P>
4 M' ]9 m& B' `$ i) O2 c+ T<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
" u) @1 J: {( c7 _: S# f<>6。龙贝格积分法</P>
<>拜托了。。。。</P>

陈了

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

西西

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

loveyuan

<>1。牛顿迭代法</P><>
  #include "stdio.h"
  #include "math.h"
0 U( X$ j3 m% |6 n1 b) p4 L  int dnewt(x,eps,js)
9 W/ ]0 ^" n: O# F& D0 c  int js;
  double *x,eps;
  { extern void dnewtf();
    int k,l;
5 H" w* b0 R' w. Y( c    double y[2],d,p,x0,x1;
6 _) x+ R; k# \8 g& ?    l=js; k=1; x0=*x;
' V/ L( I# P/ S9 I9 {9 ]    dnewtf(x0,y);
    d=eps+1.0;
$ L# T* k9 J0 t$ W/ j, ?3 `    while ((d&gt;=eps)&amp;&amp;(l!=0))
5 E7 ^1 {0 G5 H$ _5 T      { if (fabs(y[1])+1.0==1.0)
          { printf("err\n"); return(-1);}
0 G$ U4 h% F  ^# n. \        x1=x0-y[0]/y[1];
/ l  m5 k* ^1 q% C3 k+ @        dnewtf(x1,y);
        d=fabs(x1-x0); p=fabs(y[0]);
( ]0 b. ^* W# f  c2 j# @) k        if (p&gt;d) d=p;
2 i, f  t. r4 V! J        x0=x1; l=l-1;
      }
( e- |' }5 R. _4 e& p    *x=x1;
/ t/ o) v5 N. Q    k=js-l;
$ Z& Z1 s2 P9 ]' y$ G' ]    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
" e- u" Z% p3 ?4 U7 }. t2 R+ v  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
* m1 r$ ^9 K: w: F* q    js=malloc(n*sizeof(int));
  U$ n' J( i# T" ~) i1 l    for (k=0;k&lt;=n-2;k++)
      { d=0.0;
+ q+ Y/ C2 E) D8 X        for (i=k;i&lt;=n-1;i++)
# q' A" \' |6 B  t$ E( s$ n        for (j=k;j&lt;=n-1;j++)
          { u=i*n+j;
- Z; S( s6 y, z6 M7 A# l            p=ar*ar+ai*ai;
5 e  m, I1 y6 q/ E4 K* J            if (p&gt;d) {d=p;js[k]=j;is=i;}
- @7 P! q5 z4 h+ r3 s2 H6 D# f2 T6 W3 r4 w          }
- M5 x5 i6 U, A        if (d+1.0==1.0)
* f5 x/ r* O( y5 s; s          { free(js); printf("err**fail\n");
' w) U4 J: e, w            return(0);
          }
) [* l' `8 F. Z$ ?9 S, t% q+ _7 l        if (is!=k)
& {7 H8 ?( ~5 S! s' q6 }; T6 G          { for (j=k;j&lt;=n-1;j++)
* m( y' s' X, c0 F* m. X2 `& ]" e              { 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;
: b  @+ t  Q& w0 D3 S8 N' y            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
2 O, H' l) O) D/ ?; A! D          }
* W- O& l7 x  \3 N        if (js[k]!=k)
4 w' ~) q; c/ J* ?7 I5 p3 O7 b4 M          for (i=0;i&lt;=n-1;i++)
& z! ^1 _$ |% h' E; B3 y4 `1 P            { u=i*n+k; v=i*n+js[k];
              p=ar; ar=ar[v]; ar[v]=p;
- b, v/ y0 x3 Z4 x! P              p=ai; ai=ai[v]; ai[v]=p;
            }
        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
3 ?1 N: W' P& j% u8 F' O/ s4 a          { u=k*n+j;
/ T. O- d, Q6 E7 p: \5 W            p=ar*ar[v]; q=-ai*ai[v];
2 v2 Y  r- Z* a2 I; q) D0 b            s=(ar[v]-ai[v])*(ar+ai);
" c+ D3 O! x4 z            ar=(p-q)/d; ai=(s-p-q)/d;
) h  J3 }3 w: O  d( O2 F          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
  L4 W& ]% W4 |- L        s=(ar[v]-ai[v])*(br[k]+bi[k]);
) B* r$ x/ [) T0 k4 x0 X        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
* h3 _" g; _* y        for (i=k+1;i&lt;=n-1;i++)
          { u=i*n+k;
            for (j=k+1;j&lt;=n-1;j++)
* S0 m$ h. @! ^8 r0 k. j              { v=k*n+j; l=i*n+j;
                p=ar*ar[v]; q=ai*ai[v];
* d. d( _6 q2 M7 M0 l                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
4 C& [* w- w' m3 ^5 A& B, h                ai[l]=ai[l]-s+p+q;
+ v# ~1 e0 q- C/ X& V; ]              }
$ [& K- x5 H& O            p=ar*br[k]; q=ai*bi[k];
/ h' G+ J& c* o: h9 [            s=(ar+ai)*(br[k]+bi[k]);
' W% y" G4 H& ~* R& u( R* U- d8 y- B            br=br-p+q; bi=bi-s+p+q;
          }
      }
9 S/ `6 [9 P6 n    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
  m( n! d# w1 T    if (d+1.0==1.0)
& g7 [# O# f/ y) ~1 D8 h      { free(js); printf("err**fail\n");
        return(0);
      }
    p=ar*br[n-1]; q=-ai*bi[n-1];
! E# D4 F: v: y# j: P    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--)
: l/ v: V' g7 w9 ^    for (j=i+1;j&lt;=n-1;j++)
" W1 r$ W4 Z6 p9 o+ D      { u=i*n+j;
( U  @  I5 Y% \/ w        p=ar*br[j]; q=ai*bi[j];
( o$ c7 z8 X2 x        s=(ar+ai)*(br[j]+bi[j]);
        br=br-p+q;
5 E0 H9 R5 \6 s7 X2 ^        bi=bi-s+p+q;
      }
0 q8 q+ i8 G. [. _/ d* m    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
1 y+ g4 ?& x. t( c4 L# }' I      if (js[k]!=k)
        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
2 b  X) y! \& \' Y) e2 J          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
7 ^( Z0 K* L5 h: T- u9 r        }
2 H  c" S- p4 @# d  z* s2 h    free(js);
    return(1);
5 }* d/ w0 k# T7 P, P  }</P><>平方根法</P><>#include "math.h"
; Q) }6 J3 F5 f+ u5 I+ U$ B  #include "stdio.h"
5 Q/ Z# X6 n6 H% h4 w  int achol(a,n,m,d)
  int n,m;
1 b# r  X% T) ]; l% Z; E7 {  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);}
7 Z5 s. x6 v) _; u$ N# O    a[0]=sqrt(a[0]);
    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
      { u=i*n+i;
        for (j=1; j&lt;=i; j++)
) |. C8 p7 V/ g' Q4 M9 L7 G          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
8 ~3 R, e+ X$ I1 {          }
        if ((a+1.0==1.0)||(a&lt;0.0))
9 d, o3 }; O; ]6 M4 }; j6 j          { 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++)
7 o. ]5 l' a0 z9 c0 `% l6 l! V                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
/ p  h& p, O  W! |7 M% q8 P                a[v]=a[v]/a;
4 K3 U( Q) Y+ D6 @# f              }
          }
5 W' W, B, @/ s( i6 Z! V) P. A      }
    for (j=0; j&lt;=m-1; j++)
% k8 W0 x. F' ~# Q      { d[j]=d[j]/a[0];
9 R+ t+ x5 `& J; G        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
, p: q7 e1 f' P( _8 }4 Q            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;
- ^3 s4 q' ^3 |# c7 b          }
      }
    for (j=0; j&lt;=m-1; j++)
. i  W$ a( d$ e6 r% j. n      { u=(n-1)*m+j;
        d=d/a[n*n-1];
        for (k=n-1; k&gt;=1; k--)
" }/ n/ o- H4 C, }( ^) I1 l- S          { u=(k-1)*m+j;
            for (i=k; i&lt;=n-1; i++)
              { v=(k-1)*n+i;
2 k4 @3 ^5 p$ x- ]" {- [                d=d-a[v]*d[i*m+j];
              }
  x: p1 o! ]5 w1 Q) V/ a6 a            v=(k-1)*n+k-1;
  p* s, R! z0 U1 \* Q            d=d/a[v];
# T0 E# D9 U1 s( z" `          }
      }
    return(2);
& q' `& V# f5 V! t/ J1 D  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
' x- t4 S* z" ~+ h  int n;
  double x0,h,t,y[];
  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
" x# z$ U( m% O' r% T( u    z=0.0;
    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;
: p  P4 H: f" F2 B; w. D        return(z);
2 ]5 z/ `' E6 h7 z      }
    if (t&gt;x0)
( _9 T& J0 l% i$ V  E      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
      }
' N- R+ T9 t3 x$ w( T" N    else i=0;
    k=i-4;
    if (k&lt;0) k=0;
! s, p' A1 }- W    m=i+3;
: l/ a8 J: x* f4 ~, c9 E  {    if (m&gt;n-1) m=n-1;
8 g3 w3 W$ K: Z( b  b" ?    for (i=k;i&lt;=m;i++)
" o! c! G, w; a2 g# ~      { s=1.0; xi=x0+i*h;
        for (j=k; j&lt;=m; j++)
' _8 v# T6 l  a. J          if (j!=i)
            { xj=x0+j*h;
( h+ t# q6 r' ^9 J) ?6 s              s=s*(t-xj)/(xi-xj);
6 ^- z" t6 |1 ~  C            }
        z=z+s*y;
; G* H0 x9 _  H2 q$ L      }
" n  T" ]; \3 L" F( k& x    return(z);
/ ~; z+ W' u: m8 L' V# A/ |: g  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
  int n,m;
  double x[],y[],a[],dt[];
  { int i,j,k;
    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;
% d8 A2 e' N# X& W  c" t    if (m&gt;n) m=n;
    if (m&gt;20) m=20;
    z=0.0;
  K/ g! \: n) O! w* f7 q- {    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
3 ]" Y" W' k! L2 J4 ?# o# ?7 [; I    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
3 k- ~' k7 q% D1 O: ~    for (i=0; i&lt;=n-1; i++)
  E1 E) x! J  m5 ^      { p=p+(x-z); c=c+y;}
' w& ?7 W& V: h; n! H2 v    c=c/d1; p=p/d1;
- _( V: j: L( M5 y8 G    a[0]=c*b[0];
7 N" r1 [+ ~. E6 ]    if (m&gt;1)
5 h. w7 M5 s( y5 y5 ?      { t[1]=1.0; t[0]=-p;
4 ~/ J0 ?/ h0 P3 N+ n4 {; O: o7 T        d2=0.0; c=0.0; g=0.0;
# \2 @! f+ D. B        for (i=0; i&lt;=n-1; i++)
          { q=x-z-p; d2=d2+q*q;
            c=c+y*q;
6 @9 K7 I% D7 A4 a' X9 N            g=g+(x-z)*q*q;
4 s9 S  A  i- `, A$ ]( F9 W          }
  v5 ?, }, C7 E5 A4 W0 [+ c! Y0 o( z        c=c/d2; p=g/d2; q=d2/d1;
& u0 F4 Y0 e( P) ^/ `8 Z( U+ c        d1=d2;
4 Z% A& F- q9 i/ @  Q. F5 Z6 n, B        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
; N# X/ ~9 \; Q/ |      }
! Y. x# a, |  Z4 R2 _    for (j=2; j&lt;=m-1; j++)
: C' z# y. |: i8 ]1 X, o      { s[j]=t[j-1];
8 N' A5 W& a6 ]        s[j-1]=-p*t[j-1]+t[j-2];
; i2 t2 O+ X' G* z2 r' V+ n4 H7 W        if (j&gt;=3)
) F% x% k4 n- x' {          for (k=j-2; k&gt;=1; k--)
* L0 M% F0 t# l$ c; G) g- u) U- S+ V            s[k]=-p*t[k]+t[k-1]-q*b[k];
9 X7 R7 l  u7 |8 m+ a        s[0]=-p*t[0]-q*b[0];
% @, x5 _3 I$ Z$ p8 S2 I/ u6 _        d2=0.0; c=0.0; g=0.0;
5 u. w+ \( y; ]$ l" t3 G0 I        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
            for (k=j-1; k&gt;=0; k--)
              q=q*(x-z)+s[k];
            d2=d2+q*q; c=c+y*q;
            g=g+(x-z)*q*q;
+ I) F  }* c- O. a2 U          }
        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
+ l  r  y; k  N2 J0 h1 A  H        a[j]=c*s[j]; t[j]=s[j];
        for (k=j-1; k&gt;=0; k--)
          { a[k]=c*s[k]+a[k];
2 I. }  e1 C% {( g' ^            b[k]=t[k]; t[k]=s[k];
1 Q# @; e( x6 z- r" |          }
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
4 |9 E; o6 s! h7 z8 j7 ]4 ?    for (i=0; i&lt;=n-1; i++)
) ]8 w2 U# g7 E( [$ r, I9 F& F      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
  J7 V& f2 q! e          q=a[k]+q*(x-z);
        p=q-y;
        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
* a1 N8 H4 a( X! y& y9 M9 E3 X3 N4 P        dt[1]=dt[1]+fabs(p);
4 J/ T7 O7 L7 w. C/ ?      }
    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
  double fromb(a,b,eps)
  double a,b,eps;
2 ]: c. y* r1 ~1 f1 v$ }3 H. m' C  { extern double frombf();
$ L7 _: l1 V# v9 b! I7 t    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
0 a. w8 t5 j7 X" Q2 {6 p, P    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
* D: ]" Z2 j# _    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
+ Y* @7 X% K& x) [7 i) |! p7 H      { p=0.0;
        for (i=0;i&lt;=n-1;i++)
4 d# N2 W# v0 X" Q" P          { x=a+(i+0.5)*h;
4 G, G1 |' S% |5 c  u0 D            p=p+frombf(x);
( D0 b$ i) t& U8 n9 M% ^          }
& G6 P2 F3 z& F3 _6 `/ R; `7 A        p=(y[0]+h*p)/2.0;
- e$ [$ t! A) `        s=1.0;
  e' k5 W5 e! D0 U5 i# \6 v* c        for (k=1;k&lt;=m;k++)
, h/ C& K( k# r" Y  ~6 q! g5 f          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
9 K3 {/ `9 _2 \$ a- k% v            y[k-1]=p; p=q;
! O3 q$ l4 j( X% H, f          }
3 c- s0 Q8 Y5 n$ Q' N        ep=fabs(q-y[m-1]);
; O% e6 [, V7 n+ W/ Z        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
' V% U$ g2 H  w' s1 h9 M6 P6 Q" |    return(q);
  }</P><>呵呵 希望对你有用！！</P>