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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
<>1。牛顿迭代法</P>
<>1。全主消元法</P>
1 a: f) u% S) K  T" X3 P<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
% m8 m* [( y$ p<>6。龙贝格积分法</P>
" H6 ?( f( ]" M' X9 w" N: S, ]<>拜托了。。。。</P>

loveyuan

<>1。牛顿迭代法</P><>
: r& B% z0 \, [0 V$ l3 U! a( d& K  #include "stdio.h"
/ w' Q) m/ B  t- C; u  #include "math.h"
7 @4 n) I* c+ E' s  int dnewt(x,eps,js)
' s4 `: z/ y# ]9 ?) B7 n, X# r7 B9 |' e  int js;
  double *x,eps;
) W! h% T; H9 Y. c  v  { extern void dnewtf();
    int k,l;
    double y[2],d,p,x0,x1;
/ K" _. K6 }. k( [! b    l=js; k=1; x0=*x;
    dnewtf(x0,y);
    d=eps+1.0;
) t+ d- r/ ^# A7 P    while ((d&gt;=eps)&amp;&amp;(l!=0))
6 ^1 m6 Z  s6 u  X5 U7 j! l" s      { if (fabs(y[1])+1.0==1.0)
          { printf("err\n"); return(-1);}
+ {+ N% U' V. K        x1=x0-y[0]/y[1];
" w# q7 C7 u$ C9 @        dnewtf(x1,y);
7 h! p) B" U$ j0 d; b3 w        d=fabs(x1-x0); p=fabs(y[0]);
6 D9 H. L$ k! c! E+ b" ^6 x* Z        if (p&gt;d) d=p;
        x0=x1; l=l-1;
7 D3 [1 [8 l9 `0 |' l      }
    *x=x1;
* i: E* P* X1 e! ~$ J1 [) j7 O/ G    k=js-l;
/ _& }5 Z/ f# t1 \" @, P    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
* j- S5 k$ y2 ?' N  int acgas(ar,ai,n,br,bi)
9 Y6 |& n2 w; J# g: f! M% Z' n; ?' `( e  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
: p3 t$ i7 E! Y    js=malloc(n*sizeof(int));
    for (k=0;k&lt;=n-2;k++)
      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
        for (j=k;j&lt;=n-1;j++)
          { u=i*n+j;
0 h( U- G3 F1 d, s+ p( s; o) O2 `            p=ar*ar+ai*ai;
3 ~" q8 r  k2 b4 W, x            if (p&gt;d) {d=p;js[k]=j;is=i;}
6 h: Q" b6 i; \' ^          }
        if (d+1.0==1.0)
$ M( A, [# S3 R% W- t0 B9 z          { free(js); printf("err**fail\n");
- e3 I  ^+ r) z( _4 [            return(0);
2 i. F3 p4 h% l: m          }
8 n3 c- Z, M& w& i- E% X        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
                p=ar; ar=ar[v]; ar[v]=p;
& r& y2 t$ x  d9 |1 X9 Q4 A! K                p=ai; ai=ai[v]; ai[v]=p;
              }
6 m  n: L9 u+ j* A+ I2 g5 u            p=br[k]; br[k]=br[is]; br[is]=p;
+ N0 I( f' o1 h3 p5 r            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
8 u0 K/ g% x8 R2 A. n          }
        if (js[k]!=k)
+ B$ o0 @  A( ]& c6 a4 T: y- e          for (i=0;i&lt;=n-1;i++)
% L8 B3 I/ ~9 K& a0 C8 _, c( F            { u=i*n+k; v=i*n+js[k];
; Q' D+ B  X1 D3 [              p=ar; ar=ar[v]; ar[v]=p;
              p=ai; ai=ai[v]; ai[v]=p;
5 t: p% F! o) \7 g            }
( V" d! o: T- w2 Y" Y9 K4 i' }        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
1 q/ T1 _8 f' D) A          { u=k*n+j;
& c$ _8 v5 W  F, I3 K            p=ar*ar[v]; q=-ai*ai[v];
( K; J$ a8 E: ~' T( d8 H7 x7 x; v            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
4 Y0 K2 c( Y! Y, u& L9 L$ h* N          }
+ P9 e# U7 {* d' k        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
# T0 F4 [3 H2 l  _        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i&lt;=n-1;i++)
! U9 U0 ^( A1 M: {2 X/ Y$ |          { u=i*n+k;
# `! ]" M2 S! c$ D            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];
# q- l0 Y* x! E% K                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
7 n8 U6 v3 n2 w1 [: n+ A- ~7 v) X; {+ j                ai[l]=ai[l]-s+p+q;
7 c5 r2 r1 Q: S              }
2 c) n  P" B4 x            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
" z  n* d) I/ t+ `5 i' Q7 Y            br=br-p+q; bi=bi-s+p+q;
          }
      }
% e5 Q0 F6 ]6 Q  F5 s7 b& T    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
" m" M" Y0 v6 X/ z2 z+ l9 W; r6 T    if (d+1.0==1.0)
+ m3 s2 e6 f) S' B8 N      { free(js); printf("err**fail\n");
        return(0);
7 A; L5 O9 }' U. d$ ]# z* F) G, Q      }
/ L# l- {! Y5 P; o8 w6 K    p=ar*br[n-1]; q=-ai*bi[n-1];
) ~2 J* e& L3 H( M& R9 z) C4 f    s=(ar-ai)*(br[n-1]+bi[n-1]);
! J5 _, U6 ~* a: ]4 u1 ]    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
    for (i=n-2;i&gt;=0;i--)
% h1 u8 \  t" V$ q0 S    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
+ N* ?* `+ @. Y$ s        p=ar*br[j]; q=ai*bi[j];
        s=(ar+ai)*(br[j]+bi[j]);
5 p9 r; z. g6 y  n        br=br-p+q;
% V2 c5 n% \. i- e, B. E. K        bi=bi-s+p+q;
! Y/ x& m& [! i- p" p$ x& e      }
5 T. w6 X% Y- y+ z4 k    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
      if (js[k]!=k)
        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
/ V6 i/ g' p7 z2 s5 x7 c/ p, [5 A          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
5 K+ ^; O: ?" i# @        }
- ]+ Z. q1 \8 Y) z1 ?    free(js);
    return(1);
  }</P><>平方根法</P><>#include "math.h"
  #include "stdio.h"
  int achol(a,n,m,d)
- ^* w6 G% R" w7 z, [1 U9 s  int n,m;
  double a[],d[];
  { int i,j,k,u,v;
$ L& v7 g5 K: o% r- B3 p# E    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
, W7 P' |8 x/ q3 U- \      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
. {5 _2 e: G! e8 g    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
9 f4 W; I5 R# _( _      { u=i*n+i;
4 r. P9 o7 K& h$ M  ~6 n        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
. \+ {9 @0 k5 o            a=a-a[v]*a[v];
          }
        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
        a=sqrt(a);
0 G% Q5 i% b( F; X* F3 @        if (i!=(n-1))
          { 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];
4 W' O! R5 X; ]# v; S! n                a[v]=a[v]/a;
              }
          }
      }
    for (j=0; j&lt;=m-1; j++)
" f5 F# ?$ q! \4 A: B      { d[j]=d[j]/a[0];
        for (i=1; i&lt;=n-1; i++)
1 \% K( z8 y; E. u          { u=i*n+i; v=i*m+j;
4 }- L/ x% A. P1 n( b            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;
          }
      }
    for (j=0; j&lt;=m-1; j++)
+ v% G: ^8 C2 p5 G; E3 i; H      { u=(n-1)*m+j;
        d=d/a[n*n-1];
- W; b' _. a& @6 ~0 f$ C        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
. o( j0 b5 e4 [6 A3 B! @( q6 P            for (i=k; i&lt;=n-1; i++)
              { v=(k-1)*n+i;
                d=d-a[v]*d[i*m+j];
' [4 }7 _9 Z  X- z$ |5 `: N              }
            v=(k-1)*n+k-1;
            d=d/a[v];
  R' b& ^. Z' }% H3 s* N          }
2 x$ a1 \  l5 @1 @: B/ \9 D2 @. O      }
    return(2);
9 g- x* n7 h- K  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
9 Z! ~5 H5 J9 s; C4 b) B  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
    z=0.0;
    if (n&lt;1) return(z);
* h# t* |# M3 T/ w6 K( U    if (n==1) { z=y[0]; return(z);}
7 W/ E( E) X0 }" N3 [    if (n==2)
1 [' v! G  w  A7 S      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return(z);
1 W$ {: |5 \, j0 d      }
    if (t&gt;x0)
5 A# u. X  e  ^& Z5 d- p      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
      }
+ Y3 S- {4 {1 d9 L- q+ C    else i=0;
* D+ N1 J; v8 D, H    k=i-4;
    if (k&lt;0) k=0;
* V5 c  k7 G# I4 J- q/ N' H    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;
        for (j=k; j&lt;=m; j++)
6 h  _4 t: t- I4 ?7 I" Q" `          if (j!=i)
            { xj=x0+j*h;
0 Z) Z$ A* Y/ ~/ U              s=s*(t-xj)/(xi-xj);
3 U, \- d: T9 u9 L% L" Y. ^            }
        z=z+s*y;
& \7 r5 j& a, y      }
, a/ w3 r2 ?( b0 |+ T6 V( k' F$ v    return(z);
  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
" c+ f+ ?# R1 O. J* x* C  void hpir1(x,y,n,a,m,dt)
( G- W! T* [; U: s: g5 o) o  int n,m;
8 b# K1 w& A" u0 x! y  double x[],y[],a[],dt[];
  { int i,j,k;
) d; q0 q$ J# p' }( |0 n    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;
+ j. b8 u+ D6 \( q2 w    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;
    for (i=0; i&lt;=n-1; i++)
      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
4 F& F1 M8 W# ?! O8 i' c1 z    a[0]=c*b[0];
7 A! b. M7 G' z  I) ~3 Q; k    if (m&gt;1)
! A& t) Z4 ]1 W# ~: s      { t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
  }5 O* g( E0 I$ V        for (i=0; i&lt;=n-1; i++)
! T0 S- {& [3 R. A, c          { q=x-z-p; d2=d2+q*q;
            c=c+y*q;
6 Z, ?9 e3 r9 W  {            g=g+(x-z)*q*q;
1 G) l" G" I$ h          }
8 A/ \6 d# d+ Y0 ?2 ]2 {3 a        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
  d( ?5 l% j4 f9 H6 R" y        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
8 P) h3 ]9 @1 C% I+ s    for (j=2; j&lt;=m-1; j++)
      { s[j]=t[j-1];
' S* P. @& [8 W. a1 y        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
* u; Y0 s# k9 p% R4 [7 w! x          for (k=j-2; k&gt;=1; k--)
- e: e- \- K. n" 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++)
; K  ~/ s( C  I$ k! L5 ]( W          { q=s[j];
8 k  v+ t+ R8 g6 H' v* S            for (k=j-1; k&gt;=0; k--)
9 u1 U) W$ |0 K7 D              q=q*(x-z)+s[k];
            d2=d2+q*q; c=c+y*q;
            g=g+(x-z)*q*q;
          }
4 A; W5 _( y' Y6 \        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
& v; J: w* D$ t1 K, v+ G6 p: Z) Y        a[j]=c*s[j]; t[j]=s[j];
0 M4 R" j% O" L        for (k=j-1; k&gt;=0; k--)
  _5 _0 t- _& [; L' H1 R" K          { a[k]=c*s[k]+a[k];
            b[k]=t[k]; t[k]=s[k];
% P+ q- a0 q4 g# s' [' N          }
      }
  m: A& p: w! P9 u, ]    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
4 O8 z' [9 ?7 y( Y    for (i=0; i&lt;=n-1; i++)
      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
" a. m3 O: `) d) I4 q5 U2 G          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;
& X( y& c0 n. q        dt[1]=dt[1]+fabs(p);
      }
* q  G8 K% {& S6 }    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
3 |2 e* l: e( Z3 U/ l0 J8 h- g  double fromb(a,b,eps)
  double a,b,eps;
8 G' |6 r" W3 U  { extern double frombf();
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
    h=b-a;
% R4 c$ c+ P4 I4 [0 q- B    y[0]=h*(frombf(a)+frombf(b))/2.0;
. d  l( r& d4 p- C6 O- v    m=1; n=1; ep=eps+1.0;
; D  f3 |  r  I* t! g5 A    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
7 Y3 {; |: k/ S. D! E      { p=0.0;
! E1 U# d+ e, E% o! M. q- ?5 _        for (i=0;i&lt;=n-1;i++)
! v. v9 l3 V+ N          { x=a+(i+0.5)*h;
            p=p+frombf(x);
          }
        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
; L1 G0 M3 K* I3 R          { s=4.0*s;
, y! K. c% v# i0 R. V5 f, Z0 r            q=(s*p-y[k-1])/(s-1.0);
# ]  B# T# `. v4 S5 d            y[k-1]=p; p=q;
# a# _, C2 k( J, K+ m( _* @; E          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
  R3 {! T  [1 v. P' x, t    return(q);
  }</P><>呵呵 希望对你有用！！</P>

西西

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

陈了

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