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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
! [! q2 k( j! h- T<>1。牛顿迭代法</P>
! I& {# L, \/ C1 V<>1。全主消元法</P>
<>3。改进平方根法</P>
5 O$ B5 W9 L: v' \" H# _<>4。牛顿向前，向后插值法</P>
2 N4 T: C) U$ P  |" N+ }<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
7 w) D) G( Z5 N3 m: b<>拜托了。。。。</P>

loveyuan

<>1。牛顿迭代法</P><>
  #include "stdio.h"
  #include "math.h"
  int dnewt(x,eps,js)
  int js;
- z! Q2 |% S. v1 a& V: T  double *x,eps;
  { extern void dnewtf();
7 w6 ?3 T$ q$ R9 q! x( i" t    int k,l;
    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
/ N6 _( p# \4 h. a2 M# r2 R    dnewtf(x0,y);
    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);}
        x1=x0-y[0]/y[1];
        dnewtf(x1,y);
- |3 q# y% o. Z$ p7 A8 b: V        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
      }
# H% a  c  a! p) A    *x=x1;
5 O$ z  f$ A* V4 E% Q    k=js-l;
    return(k);
  T$ {- N9 N, N4 N0 L  }</P><>全主消元法</P><>#include "stdlib.h"
" p7 y% B5 `2 `2 G# h; y" ~6 T, u  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
% v4 f! L5 O! G& k  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
! m4 m8 C' I- }2 F$ p- R; g    double p,q,s,d;
; N9 B& P/ q6 I; m5 p    js=malloc(n*sizeof(int));
' t2 ~, I: R) F; i) {/ `    for (k=0;k&lt;=n-2;k++)
* L8 U. S1 h) v- C4 }- e      { d=0.0;
2 r0 q. M# ]6 {; _; a8 @: Y. D        for (i=k;i&lt;=n-1;i++)
. e" e8 K/ k! m; z1 o        for (j=k;j&lt;=n-1;j++)
          { u=i*n+j;
/ {! Q" x! p0 S# Z& W( `            p=ar*ar+ai*ai;
            if (p&gt;d) {d=p;js[k]=j;is=i;}
; S6 V. }6 [  l) Z' O          }
# _+ K- H( M5 j; f3 B* @        if (d+1.0==1.0)
' i7 }, `' K* B6 ~5 [          { free(js); printf("err**fail\n");
# b/ V: E  U3 w6 S, j! M( M) y            return(0);
          }
        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
# D/ C6 P* `- b/ h/ _                p=ar; ar=ar[v]; ar[v]=p;
                p=ai; ai=ai[v]; ai[v]=p;
6 B9 H6 _- |; k              }
            p=br[k]; br[k]=br[is]; br[is]=p;
: h8 T8 [  c3 K            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
0 k; m& g7 h% Z8 R$ @          }
        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
' }! ]( Z8 S- u( ]. N, |( x6 V            { u=i*n+k; v=i*n+js[k];
) Q) ~' j; L! E, j) H+ K+ o, R+ I# f              p=ar; ar=ar[v]; ar[v]=p;
              p=ai; ai=ai[v]; ai[v]=p;
            }
        v=k*n+k;
$ e! u" v$ J: w# l$ o        for (j=k+1;j&lt;=n-1;j++)
5 O& c% O/ n, l          { u=k*n+j;
            p=ar*ar[v]; q=-ai*ai[v];
6 k% {' ~, \0 y. |( M            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
, m- T& Z' {) i  l0 o6 t$ d# z          }
, ]# h1 \6 Z& O+ 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;
1 R2 [6 k: M% d        for (i=k+1;i&lt;=n-1;i++)
          { u=i*n+k;
0 k$ f2 q7 m) Q# Y6 S            for (j=k+1;j&lt;=n-1;j++)
              { v=k*n+j; l=i*n+j;
                p=ar*ar[v]; q=ai*ai[v];
                s=(ar+ai)*(ar[v]+ai[v]);
/ i# U9 }0 @2 i( Z8 O* n                ar[l]=ar[l]-p+q;
, x- p5 Z9 S. \% \, v                ai[l]=ai[l]-s+p+q;
8 }. R: F) g/ K& G" T* }0 l              }
+ Y1 q6 z/ O; e            p=ar*br[k]; q=ai*bi[k];
; B- p, U9 j& b            s=(ar+ai)*(br[k]+bi[k]);
  H4 A" j" P5 e* B            br=br-p+q; bi=bi-s+p+q;
% X& I+ r2 @! E; l+ d4 X          }
      }
' b/ d2 e/ ]) n: v8 @0 S) D( z8 r    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
        return(0);
; F% o. L2 p  \; [$ e: [      }
. G& T% |9 {- }9 n( H    p=ar*br[n-1]; q=-ai*bi[n-1];
9 d6 `8 P! c6 ]    s=(ar-ai)*(br[n-1]+bi[n-1]);
; t, j; Y. G/ D+ d, x. ^% V    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
" S$ j) Y! K* }    for (i=n-2;i&gt;=0;i--)
. k8 ~1 F) t$ J% k# Z    for (j=i+1;j&lt;=n-1;j++)
1 a, P0 y  z5 H+ i      { u=i*n+j;
        p=ar*br[j]; q=ai*bi[j];
        s=(ar+ai)*(br[j]+bi[j]);
; |' |4 C  w" W" r3 U/ y        br=br-p+q;
+ u- Z- Y+ X4 d5 V        bi=bi-s+p+q;
      }
0 l; r& l- c0 |6 F    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
7 B' r, N+ q: I3 [      if (js[k]!=k)
, I2 b5 [9 k/ T) b8 p0 I6 j        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
; b4 F5 H. H5 L. i& |          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
        }
    free(js);
3 ^% w9 E8 I9 U5 C# G    return(1);
  }</P><>平方根法</P><>#include "math.h"
4 n, C% U( d6 X5 [- c# e  #include "stdio.h"
  int achol(a,n,m,d)
  int n,m;
# \" v- l* q* O  double a[],d[];
; [: P5 R! ~- n  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
      { printf("fail\n"); return(-2);}
, Q" u  Y. t( X, R  ^! a6 S% t    a[0]=sqrt(a[0]);
  U* u7 q5 u- T1 N3 d4 I    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
+ l8 i& l" y2 D      { u=i*n+i;
        for (j=1; j&lt;=i; j++)
- P- P4 j& h- H6 G" E2 @5 Q          { v=(j-1)*n+i;
1 |& j8 g+ u* B; \6 I# O3 @            a=a-a[v]*a[v];
          }
        if ((a+1.0==1.0)||(a&lt;0.0))
5 I* B. Z" u" R5 A0 o          { printf("fail\n"); return(-2);}
: v# c0 ?2 X7 j$ ?$ [8 m5 y        a=sqrt(a);
        if (i!=(n-1))
; F* K/ _1 r% B: i0 w* o4 J          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
                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;
$ t2 D( S2 I% n) ^8 M- m/ A; q              }
          }
      }
    for (j=0; j&lt;=m-1; j++)
      { d[j]=d[j]/a[0];
2 }  R, a" R9 P2 k1 k        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
7 }; x5 V( H, v  ]9 J' k            for (k=1; k&lt;=i; k++)
3 |7 f) o) n! ]: p. e              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a;
+ ?" @& ?  I" q) ~4 U$ l          }
      }
    for (j=0; j&lt;=m-1; j++)
9 \5 N, u- l; Y      { u=(n-1)*m+j;
( \9 v3 Z, S: }        d=d/a[n*n-1];
        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
6 k9 E( B0 |! u& ~  H8 ~) N1 f6 l            for (i=k; i&lt;=n-1; i++)
" d6 _% [2 d! q1 m; j3 V4 j              { v=(k-1)*n+i;
& y5 H" m# t6 C- F, Q' l2 H                d=d-a[v]*d[i*m+j];
              }
            v=(k-1)*n+k-1;
2 k( `# |# j9 l' O' L% C            d=d/a[v];
' s0 k( `1 d0 C' d$ Q          }
      }
    return(2);
  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
. e+ t$ {0 j! ^- L% _+ h# g4 q  double x0,h,t,y[];
  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
    z=0.0;
+ f) x- O/ o) |/ ~/ }' @    if (n&lt;1) return(z);
    if (n==1) { z=y[0]; return(z);}
/ a& ^; ]& z, G$ o1 S* q9 X$ t    if (n==2)
* j" A* z6 z% h" }      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
  g* i& }* Q0 F) \# X" F1 O( D7 v; H' W        return(z);
      }
    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
& B2 [) H- S8 s& i      }
: E$ w  P/ V2 k0 z& o* ^! z/ ]" I    else i=0;
    k=i-4;
    if (k&lt;0) k=0;
( }4 L# [3 f3 [7 l& P/ P% l    m=i+3;
, A8 {  Z) B# N* w& \8 n    if (m&gt;n-1) m=n-1;
    for (i=k;i&lt;=m;i++)
      { s=1.0; xi=x0+i*h;
        for (j=k; j&lt;=m; j++)
2 B' @4 F6 {2 Z/ d. I          if (j!=i)
/ u9 g$ d% N( T            { xj=x0+j*h;
  H9 K" F, j8 C  g% N+ [              s=s*(t-xj)/(xi-xj);
! ^0 w' u* Y6 i            }
        z=z+s*y;
      }
0 L7 y% B  ]% x+ v) n    return(z);
" P6 J9 [2 o, _* I& Q  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
  int n,m;
6 c! w) W+ k6 j  double x[],y[],a[],dt[];
  { int i,j,k;
6 k( b8 o: K7 \: J    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;
9 l; b8 q( ^& l- n    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;
* r; F, b4 `; S4 L* ~( ^' U8 M% ]    for (i=0; i&lt;=n-1; i++)
      { p=p+(x-z); c=c+y;}
8 u8 W: i9 A/ e* P9 {    c=c/d1; p=p/d1;
5 e* D& j* J: ?3 P0 C    a[0]=c*b[0];
- l* [9 u( p3 F5 b/ R- d+ o    if (m&gt;1)
      { t[1]=1.0; t[0]=-p;
7 W$ j. K) T( O6 ?0 D9 B        d2=0.0; c=0.0; g=0.0;
. e) n3 w' }0 u4 c0 h        for (i=0; i&lt;=n-1; i++)
8 s% I/ f3 n1 Q2 ?          { q=x-z-p; d2=d2+q*q;
2 d2 Y* d0 s) o( x            c=c+y*q;
+ a" k' S& F  |) g+ K2 g            g=g+(x-z)*q*q;
          }
        c=c/d2; p=g/d2; q=d2/d1;
: h  t& T2 B+ ]! ?8 ]        d1=d2;
        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
" @' Q& j6 t$ |2 @7 Z    for (j=2; j&lt;=m-1; j++)
" X: y3 ^4 a0 u, {) _. e      { s[j]=t[j-1];
9 o% {8 f  k4 I        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
1 h) O/ `$ ^. h: o          for (k=j-2; k&gt;=1; k--)
! X. R: @1 l5 j5 b            s[k]=-p*t[k]+t[k-1]-q*b[k];
        s[0]=-p*t[0]-q*b[0];
        d2=0.0; c=0.0; g=0.0;
        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;
3 B! q4 F% f: ?2 a            g=g+(x-z)*q*q;
          }
        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];
& [1 ?8 f4 n+ x8 m0 m, i/ p            b[k]=t[k]; t[k]=s[k];
4 t+ b' |$ d5 R3 `. u6 J- L          }
      }
* T' e$ f* W, V; X. x( z8 l    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
    for (i=0; i&lt;=n-1; i++)
      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
          q=a[k]+q*(x-z);
        p=q-y;
2 q' G* D+ n$ y9 `3 x        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
9 W; g. C3 T1 H$ u        dt[1]=dt[1]+fabs(p);
" O! `  f* o5 X0 n- d4 m9 T      }
    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
; d8 A6 f1 ~2 |( D; ~0 \: R  double fromb(a,b,eps)
  double a,b,eps;
! C' [/ J& P8 D: [  { extern double frombf();
- {, Y. l8 y- |' R4 a    int m,n,i,k;
" J/ x( ]$ }3 h3 V2 G5 g  n9 k    double y[10],h,ep,p,x,s,q;
    h=b-a;
) I6 S; \0 s) P9 R    y[0]=h*(frombf(a)+frombf(b))/2.0;
8 J$ D' O0 b. Y8 ~# P5 K( h3 [7 n    m=1; n=1; ep=eps+1.0;
8 a* N& u% J& K- S    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
3 b; X( Q0 M6 ^* V- q      { p=0.0;
# l$ l: R6 A/ l' {% _. T        for (i=0;i&lt;=n-1;i++)
1 p' A% _9 c+ O) b( _8 H; q          { x=a+(i+0.5)*h;
            p=p+frombf(x);
/ m. L& g. _4 y" d4 p2 m9 ?          }
        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
; P) C$ G/ L( x2 S* D) {) d. m; K          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
8 ], \! L0 f0 y& N6 f0 @            y[k-1]=p; p=q;
" y# Y0 a6 D: |( l; L          }
! n0 ^2 {8 u: x        ep=fabs(q-y[m-1]);
6 _5 L7 |) e+ a0 G        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
    return(q);
  }</P><>呵呵 希望对你有用！！</P>

西西

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

陈了

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