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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
5 s' c% @* u7 I4 n5 ?/ i<>1。牛顿迭代法</P>
8 N: i4 a! g* @& z$ w" P. R<>1。全主消元法</P>
<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
6 E8 ~5 {8 c. h. n$ l) D, c, Q<>6。龙贝格积分法</P>
8 q' V# ~# T2 h2 y- R; u  d<>拜托了。。。。</P>

loveyuan

<>1。牛顿迭代法</P><>
9 |& e% I2 z7 i8 P4 r  #include "stdio.h"
  #include "math.h"
  int dnewt(x,eps,js)
  Q9 P6 T4 a1 a/ ]  int js;
  double *x,eps;
  { extern void dnewtf();
: z9 ~( F9 W. `- l% e    int k,l;
  H2 _# g0 o. ]    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
- X" |1 `5 r( L4 j( a) I$ K    dnewtf(x0,y);
    d=eps+1.0;
; X) @, Y& a# x, h- M4 G    while ((d&gt;=eps)&amp;&amp;(l!=0))
$ E) |* f# o/ B' Z1 U      { if (fabs(y[1])+1.0==1.0)
4 j$ n8 {* U+ q          { printf("err\n"); return(-1);}
        x1=x0-y[0]/y[1];
0 a! [2 j- e- J        dnewtf(x1,y);
        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
9 L0 Z+ o: {/ p0 w% e5 y      }
    *x=x1;
    k=js-l;
    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
8 I4 o# H4 a/ g  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
( `) n! z6 t- w$ a; X& l0 j2 Q    double p,q,s,d;
    js=malloc(n*sizeof(int));
. r6 d( D- ], J/ a& O, f& {/ w    for (k=0;k&lt;=n-2;k++)
, m5 ^: l) o6 L, g2 D      { d=0.0;
- g0 ~0 ?' R! Y( c        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;}
          }
  i9 g# b+ P! {% \        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
# j3 u- k! p. V            return(0);
7 B: N6 L8 }/ G  l: N1 _          }
        if (is!=k)
3 O- y7 l, j7 a  O% P1 s; y          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
                p=ar; ar=ar[v]; ar[v]=p;
: ~9 D) T, v) Q* T3 K; J                p=ai; ai=ai[v]; ai[v]=p;
              }
/ N9 _$ U" B% ^' ^* Q, h7 c            p=br[k]; br[k]=br[is]; br[is]=p;
            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
- T' s% z6 x: ~# O+ r          }
        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
/ R+ y0 A3 Q5 C            { u=i*n+k; v=i*n+js[k];
; d8 S$ V. s9 n( R' D" P' F              p=ar; ar=ar[v]; ar[v]=p;
* E5 I, f$ {/ S+ G6 N" F* b              p=ai; ai=ai[v]; ai[v]=p;
% N/ l1 [. h2 }3 |* c            }
        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
* E; w+ \# v, f$ @8 j; t$ _( M          { u=k*n+j;
. t+ ^( V7 V* J; ]8 L4 C# m2 M& r% A            p=ar*ar[v]; q=-ai*ai[v];
            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
- q  J9 A1 h- i4 u4 G          }
        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;
6 m& ^7 {1 A6 V        for (i=k+1;i&lt;=n-1;i++)
( h7 |. x( X  A2 ?          { u=i*n+k;
& w7 P2 ~3 {2 E            for (j=k+1;j&lt;=n-1;j++)
+ L: Y& [4 J/ ^7 Q) N( R              { v=k*n+j; l=i*n+j;
                p=ar*ar[v]; q=ai*ai[v];
2 o4 B- o+ n3 Y5 E                s=(ar+ai)*(ar[v]+ai[v]);
% ?' E7 S7 F* L1 w  x                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;
          }
      }
    u=(n-1)*n+n-1;
9 X8 j5 o( c/ L* h2 \    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
- T9 X$ e; T% {0 t        return(0);
      }
; \( W$ k- Q0 B# `    p=ar*br[n-1]; q=-ai*bi[n-1];
    s=(ar-ai)*(br[n-1]+bi[n-1]);
  `; j, J2 Y) d2 V    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
, S9 N  T3 Y0 M( t    for (i=n-2;i&gt;=0;i--)
    for (j=i+1;j&lt;=n-1;j++)
9 q9 x" B* {4 |* z! G! H      { u=i*n+j;
        p=ar*br[j]; q=ai*bi[j];
9 L: ]% J; t6 G, \9 }0 l        s=(ar+ai)*(br[j]+bi[j]);
8 J8 c# t5 x+ W" M        br=br-p+q;
        bi=bi-s+p+q;
      }
    js[n-1]=n-1;
! \# ]" o. F4 Q9 b5 r1 U. G# i    for (k=n-1;k&gt;=0;k--)
      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;
2 s( B; x. V" j6 s7 \# V' I; A6 q% v9 n* J        }
  o, \( V% w; I* r, S/ k) _8 [    free(js);
    return(1);
8 q  Y% n2 [% f+ U# y( w+ D3 @0 b  }</P><>平方根法</P><>#include "math.h"
  #include "stdio.h"
2 W" e3 d3 V. k6 d0 B* D- ^. E' ]" q7 l  int achol(a,n,m,d)
) E, K5 M% ?' Q  int n,m;
  double a[],d[];
, P) T* q  b: ]# f& c: b3 p0 ^3 S  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
# K. I( J, i. y6 G% v  m7 n      { printf("fail\n"); return(-2);}
" F8 W# U, J4 W  G9 ]    a[0]=sqrt(a[0]);
. Q! e5 W+ }) M6 N! \% `    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
6 E2 A. X+ j8 K      { u=i*n+i;
        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
          }
        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
        a=sqrt(a);
% c9 V5 ?& k# {/ h4 X6 Y        if (i!=(n-1))
+ i9 u3 Z: v+ e5 d+ f3 @* |          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
$ p2 N4 T( e: K8 s" d& d                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;
! ^$ Y/ c/ |- `( O4 g( E              }
          }
      }
8 N$ t' U* y% l2 Y/ r5 A    for (j=0; j&lt;=m-1; j++)
+ J, Y6 X! I  B4 I      { d[j]=d[j]/a[0];
$ b) `' y8 a: ]! n        for (i=1; i&lt;=n-1; i++)
" C. E& y# b, m          { u=i*n+i; v=i*m+j;
            for (k=1; k&lt;=i; k++)
. j% h+ m0 ?+ [              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a;
          }
; y8 M: p0 q  ]3 s! R  O      }
) F5 l# z1 j0 A; U) W. r( q: n' D    for (j=0; j&lt;=m-1; j++)
' T; Q7 b- H- O      { u=(n-1)*m+j;
' @, b2 e- c2 ?! S$ H1 ^        d=d/a[n*n-1];
& C* h5 I( M1 |$ w# x1 A/ J  j        for (k=n-1; k&gt;=1; k--)
/ v. @* R7 Z+ f% P) a, e0 Z          { u=(k-1)*m+j;
, K2 h) [; y, p( L  @* t$ P* {3 V            for (i=k; i&lt;=n-1; i++)
  ^; ^& f' Y+ c              { v=(k-1)*n+i;
                d=d-a[v]*d[i*m+j];
% z- P  }9 h, d, h0 P' @' v              }
: b8 b. _+ \1 ~( \            v=(k-1)*n+k-1;
            d=d/a[v];
          }
+ I, L( h! V! [, z      }
    return(2);
  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
  { int i,j,k,m;
; _; l# g2 b) e3 y/ C# P4 S2 w6 N    double z,s,xi,xj;
, x) k5 w9 t3 f/ n6 W) X+ R    float p,q;
    z=0.0;
# S0 }$ a1 R1 ]0 R    if (n&lt;1) return(z);
    if (n==1) { z=y[0]; return(z);}
9 w0 ~  a+ S) W3 q6 B4 l2 J7 Y    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
. p4 D4 E; z+ d' E8 U( [( ?        return(z);
      }
3 r) w7 |4 d5 |3 }    if (t&gt;x0)
( K+ I' D! f: k: e# N      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
6 e( D& \/ j% X% V, B      }
' \' _* C( d# E" {$ k    else i=0;
    k=i-4;
    if (k&lt;0) k=0;
0 S& ^* j! \( x    m=i+3;
    if (m&gt;n-1) m=n-1;
- p4 J* D$ z# O, Q1 V' s    for (i=k;i&lt;=m;i++)
      { s=1.0; xi=x0+i*h;
        for (j=k; j&lt;=m; j++)
- M( O* s$ \! P$ `7 U          if (j!=i)
4 J* ~9 `% B7 K/ j$ d  i* w            { xj=x0+j*h;
6 b% r7 ]5 d  x* O/ c0 ]              s=s*(t-xj)/(xi-xj);
/ W( H1 k9 y8 _2 {" e1 @2 `! W            }
% F- L: A! d2 S! F        z=z+s*y;
      }
/ ]5 c; j1 V; M/ \# y' u4 C  C    return(z);
  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
, r$ t; o+ J% H% p# P' r! N: l  void hpir1(x,y,n,a,m,dt)
* o" h3 g# B- q2 \7 u  int n,m;
  double x[],y[],a[],dt[];
% y: j/ [5 M/ _6 ?  { int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
! n6 Y+ B0 ^! j! d$ ?    for (i=0; i&lt;=m-1; i++) a=0.0;
    if (m&gt;n) m=n;
    if (m&gt;20) m=20;
5 F- ^8 N2 a7 R# d    z=0.0;
    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
4 h9 x; v2 Y) N8 ^+ v# j' o    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
    for (i=0; i&lt;=n-1; i++)
% @, ]! \" J- S+ v      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
. @0 `5 ^% h# j; n- w9 O) |* M2 l    a[0]=c*b[0];
) `5 R7 l; g$ B. T2 d+ T    if (m&gt;1)
+ p) d6 O3 @; I- h. N      { t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
( q* P  d+ F; v% c: H5 Q" U* V/ ~        for (i=0; i&lt;=n-1; i++)
          { q=x-z-p; d2=d2+q*q;
% G" M* q6 t- t2 ~- N9 Y8 G2 ~5 Q            c=c+y*q;
            g=g+(x-z)*q*q;
0 t5 [7 |% h8 P( K! z          }
        c=c/d2; p=g/d2; q=d2/d1;
* z" C- u! ]+ O9 Y( E        d1=d2;
        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
7 q$ F+ `7 F9 ?5 z- m4 l$ x4 N      }
% r! G0 j$ n; B* l. {  {    for (j=2; j&lt;=m-1; j++)
# X  r2 Z( _9 i+ |      { s[j]=t[j-1];
! s. C) S. Q! Q8 x        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
/ a0 _# ~0 ?$ ?4 n& w          for (k=j-2; k&gt;=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
5 {6 y7 y5 C( p        s[0]=-p*t[0]-q*b[0];
        d2=0.0; c=0.0; g=0.0;
* r/ k. E, x8 M% A; N        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
1 F% {: o: R; Y* q. A6 [            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;
          }
; |/ Y" M4 k! \- K) p  D$ Y        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
        a[j]=c*s[j]; t[j]=s[j];
        for (k=j-1; k&gt;=0; k--)
          { a[k]=c*s[k]+a[k];
3 a7 h$ o) B& m            b[k]=t[k]; t[k]=s[k];
/ }) H% O4 r8 p) h          }
$ u- f. B* {2 \9 `8 [0 z# I! _6 n      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
( d9 j1 ?* x! J( M; P1 p- d    for (i=0; i&lt;=n-1; i++)
! a! _# k% f, b& N      { q=a[m-1];
) N3 G1 a, \- C2 ~* g- `% F: @' Z        for (k=m-2; k&gt;=0; k--)
          q=a[k]+q*(x-z);
0 ^3 V* H, ?, ^2 l) l        p=q-y;
+ ]1 e- _# O( F; M) q        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
4 G, t- N6 H# T- a2 J" t. }& E        dt[0]=dt[0]+p*p;
7 z  {1 J$ o9 i6 o7 h1 d: [7 ?        dt[1]=dt[1]+fabs(p);
7 N* Q  a2 s; ?9 L5 R: |      }
    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
/ v5 v, e. z, B" w  double fromb(a,b,eps)
  double a,b,eps;
  { extern double frombf();
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
0 G3 ^+ @+ j# }& `* X7 {  V; ?    h=b-a;
  {# ^% r+ f# z2 l9 e  d$ \5 d    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
+ `6 S7 C* q) l1 ^5 s3 B    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
      { p=0.0;
& b2 D- s: w7 d        for (i=0;i&lt;=n-1;i++)
3 z- h! I" O5 |" J8 P          { x=a+(i+0.5)*h;
            p=p+frombf(x);
          }
' L$ x' U3 u, O' L        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
, Q. N0 f7 y: p$ `& F3 ^          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
. r0 U( S; P2 ]; q! U            y[k-1]=p; p=q;
* l. D. ^0 T; @          }
: j9 w2 I4 F2 s5 d; \        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
$ n' m0 f+ h" o" r    return(q);
  }</P><>呵呵 希望对你有用！！</P>

西西

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

陈了

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