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

西西

<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
4 R' {" w" w8 L* ^% \9 {<>1。牛顿迭代法</P>
<>1。全主消元法</P>
<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
9 d: @% @: a, |" y6 |& ?3 o<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
5 `  R9 s: n$ |% a4 i1 l. X<>拜托了。。。。</P>

loveyuan

<>1。牛顿迭代法</P><>
' j0 f5 o2 m: ^- C6 `: Q3 L3 B  #include "stdio.h"
  #include "math.h"
6 \7 s1 @3 ^) ~  int dnewt(x,eps,js)
  int js;
  M8 z2 H5 M, w. a! l; c  double *x,eps;
  { extern void dnewtf();
1 K8 o6 a/ }1 K    int k,l;
$ H% s& X. F5 t    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
3 S8 M% p4 |- V+ [1 W    dnewtf(x0,y);
  c2 x% P1 G$ v; @9 z    d=eps+1.0;
% ?; p* r4 F# {* ?1 q    while ((d&gt;=eps)&amp;&amp;(l!=0))
" s1 P' r) H1 \2 D% |% k# [      { if (fabs(y[1])+1.0==1.0)
' o* N* f/ \: P) e" ]( g" d- ]          { printf("err\n"); return(-1);}
        x1=x0-y[0]/y[1];
; f1 D2 q& P) _* X/ |5 X6 c        dnewtf(x1,y);
2 o- v5 ^- Y9 G( T+ j        d=fabs(x1-x0); p=fabs(y[0]);
9 f+ `/ O) }3 ^) T        if (p&gt;d) d=p;
; M! h+ w) b% t        x0=x1; l=l-1;
      }
    *x=x1;
1 N$ M% q' @, m3 i% H& b    k=js-l;
+ l' L# V9 ^; |8 D! T  I    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
) x0 O% K6 @# \  #include "stdio.h"
7 s! M% z* v- j* u  int acgas(ar,ai,n,br,bi)
2 N' p  c, G& r8 {% U: |. Y  int n;
  double ar[],ai[],br[],bi[];
( {. B0 k( B6 [2 {* J+ z/ }  { int *js,l,k,i,j,is,u,v;
1 ]; l3 [: x) E* f' d9 z    double p,q,s,d;
    js=malloc(n*sizeof(int));
+ p2 r( m; e* P" `1 F    for (k=0;k&lt;=n-2;k++)
      { d=0.0;
# _* w3 Z  R! a+ G) u5 m        for (i=k;i&lt;=n-1;i++)
        for (j=k;j&lt;=n-1;j++)
! _9 {1 I& x( |          { u=i*n+j;
: c5 ]& A+ d+ V0 Z3 v# N            p=ar*ar+ai*ai;
' F' L" N! @2 d- j7 R8 `            if (p&gt;d) {d=p;js[k]=j;is=i;}
) i* @) u) [0 p% Y          }
) }; H$ {( |$ k6 ?9 E8 a3 D        if (d+1.0==1.0)
& t: B0 D1 Q8 P$ o          { free(js); printf("err**fail\n");
/ s; k- n- }( `9 K2 a9 f3 u: Y            return(0);
          }
  r: }/ W# `+ m        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
( b8 ~5 @) V  r- q0 k& A              { u=k*n+j; v=is*n+j;
# z9 ~7 g* ^) D7 S, _- L                p=ar; ar=ar[v]; ar[v]=p;
5 X' R$ Y0 h; S1 g                p=ai; ai=ai[v]; ai[v]=p;
( d1 P: M0 Q# E1 I# G  v              }
            p=br[k]; br[k]=br[is]; br[is]=p;
+ x1 n2 X0 m( k3 a* q( b            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
# o. w' Y" G' S, j# P5 q+ ^4 L          }
1 r# M' |/ m& l        if (js[k]!=k)
/ N7 s' a* y5 s+ t# ]$ X4 d! s          for (i=0;i&lt;=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
              p=ar; ar=ar[v]; ar[v]=p;
              p=ai; ai=ai[v]; ai[v]=p;
            }
* x( C5 ~# N. l' y        v=k*n+k;
9 P. M$ ^+ P& v; O  I! D% l4 H+ S7 G        for (j=k+1;j&lt;=n-1;j++)
          { u=k*n+j;
; V$ h7 g0 d9 y3 }& k9 t- V# a: H            p=ar*ar[v]; q=-ai*ai[v];
# v# ^! v7 F) d% v# u: Q3 ?) K$ E2 H            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
; n" E* E1 W( [2 R' A          }
& ]) \0 m0 L$ g  }9 C0 I        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;
        for (i=k+1;i&lt;=n-1;i++)
% \! }7 G- I1 @/ d          { u=i*n+k;
            for (j=k+1;j&lt;=n-1;j++)
  M+ \  c7 v0 O0 j/ v! U- T- Y              { v=k*n+j; l=i*n+j;
0 J# y# B. e9 X1 a                p=ar*ar[v]; q=ai*ai[v];
5 l: R1 D" z- S/ d% ^8 ]% B* z                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
            p=ar*br[k]; q=ai*bi[k];
+ `+ q# X; k2 \% k8 t- r            s=(ar+ai)*(br[k]+bi[k]);
            br=br-p+q; bi=bi-s+p+q;
+ R" E2 _, B! I1 N; A. B! A          }
      }
( ^. K4 J! ~3 W7 K9 `    u=(n-1)*n+n-1;
3 I3 ^# M, [7 L; K, s    d=ar*ar+ai*ai;
% k* g; i, W! u1 i( j    if (d+1.0==1.0)
: j$ e9 q+ u0 }  X" u      { free(js); printf("err**fail\n");
' E& U- `+ f" K& T/ o        return(0);
+ q0 o7 F3 {4 U0 B      }
    p=ar*br[n-1]; q=-ai*bi[n-1];
) x7 [& C' ?, [2 h' `% e    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--)
. @, w' d' V' U$ A0 P3 ]% q( K    for (j=i+1;j&lt;=n-1;j++)
" X' m9 k* p* D5 }$ G1 O6 n8 I      { u=i*n+j;
# R* \" h( O; V# t& Z" E2 G4 B0 a        p=ar*br[j]; q=ai*bi[j];
        s=(ar+ai)*(br[j]+bi[j]);
- \2 w' G- j) W% G        br=br-p+q;
        bi=bi-s+p+q;
! h; ]2 q! f- Z+ s: Z  _0 V      }
: C/ t# V5 D- Y/ ], {" f    js[n-1]=n-1;
    for (k=n-1;k&gt;=0;k--)
      if (js[k]!=k)
6 @# L& t/ O1 U. |3 _4 {        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
* x3 o! q3 R2 i- x( V8 Q3 C- M7 a        }
2 S/ x# H; z2 C) S4 L# y    free(js);
+ d0 Z/ h  F, e- O) r    return(1);
  }</P><>平方根法</P><>#include "math.h"
2 {( Y$ V" }# t* I  P  #include "stdio.h"
# A% {& ?( D+ Q  u1 V3 m& E! j4 d  int achol(a,n,m,d)
  int n,m;
- d" O4 s1 ?& f) n4 d% j  double a[],d[];
  W; o! h! q& k0 F6 s' B# i/ b  E  { int i,j,k,u,v;
6 I" {; M! W& G: d    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
* O* U9 {# I* v0 `" @3 n    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
2 _/ S. n6 l' U' X8 {    for (i=1; i&lt;=n-1; i++)
      { u=i*n+i;
4 b5 y+ J' f5 Q        for (j=1; j&lt;=i; j++)
" s: ?# H8 x+ J  _0 E8 _( L- @          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
          }
        if ((a+1.0==1.0)||(a&lt;0.0))
! h1 y- ~9 L, L4 ]' t          { printf("fail\n"); return(-2);}
        a=sqrt(a);
- |% S7 h/ P& T! V; B6 ]        if (i!=(n-1))
+ X5 A. q* O) c+ J6 N! T          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
7 g# D+ n+ F0 n+ ^% q                for (k=1; k&lt;=i; k++)
% e* ?4 z0 x7 ]% t5 Q3 x7 J; k* ]                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
                a[v]=a[v]/a;
- u" T2 G6 y" v# h) K              }
3 r% [) ]* X- Z5 n* P! k          }
. }5 |6 q$ l+ Q- q7 d* k( n      }
+ D1 h' K. e; c( c) M3 w( p    for (j=0; j&lt;=m-1; j++)
      { d[j]=d[j]/a[0];
        for (i=1; i&lt;=n-1; i++)
5 [1 v8 {& w' Q2 f: u1 U  |# o          { u=i*n+i; v=i*m+j;
; I: E* z) I) W" C. l7 q            for (k=1; k&lt;=i; k++)
9 w7 \1 w8 t& l: f& G. M8 X8 O              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
% P: n- G" S0 [            d[v]=d[v]/a;
4 n/ {6 A! H4 j. F6 ~' s          }
' Z' I# v! T- S  ~. c9 V- q- b      }
+ l: q1 l9 [  W; i5 ^; e3 j1 V    for (j=0; j&lt;=m-1; j++)
6 ?* H" v) J6 v- O      { u=(n-1)*m+j;
        d=d/a[n*n-1];
        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
            for (i=k; i&lt;=n-1; i++)
( G5 h3 d9 V1 Z0 x+ X" k( Z              { v=(k-1)*n+i;
                d=d-a[v]*d[i*m+j];
* g% m9 \: ^* Q' O$ k5 l; [              }
* D, e! r% Z+ G/ @5 o8 q2 y            v=(k-1)*n+k-1;
5 l. D: E. i& h# |; s            d=d/a[v];
7 r" ^9 T. a( `; c. e7 ~, b          }
9 {( ?1 I6 H* ]      }
    return(2);
- m: {4 S6 s& D  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
    z=0.0;
    if (n&lt;1) return(z);
    if (n==1) { z=y[0]; return(z);}
7 ]$ B. s1 M8 I2 j/ t! M    if (n==2)
$ b2 g7 w- a& u7 t. P      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
) t) |" I( ]' }        return(z);
5 `3 K, {' f- A1 @% E" @      }
    if (t&gt;x0)
4 k# ~) g# l) c. e" K      { p=(t-x0)/h; i=(int)p; q=(float)i;
  W& U, Q5 p$ _7 Y5 a5 ]6 G        if (p&gt;q) i=i+1;
. G4 D0 E% m3 z) h7 E- D( L, U/ U      }
    else i=0;
' q: _, Q- B( T2 n3 V3 o    k=i-4;
% ]/ _& |! s6 _& P1 r. |    if (k&lt;0) k=0;
    m=i+3;
  T8 L; K6 c: }- ^. F# o$ V    if (m&gt;n-1) m=n-1;
    for (i=k;i&lt;=m;i++)
5 y# X: r% i: d; O      { s=1.0; xi=x0+i*h;
        for (j=k; j&lt;=m; j++)
# N  S# n1 y* {5 }9 ~0 x          if (j!=i)
$ [! s( i/ d: y# O            { xj=x0+j*h;
5 D( i. [0 [$ q) F- ~8 f* g              s=s*(t-xj)/(xi-xj);
) m1 n! W" @7 }. L7 @' T            }
8 c8 q3 Z( P* C) `+ u% f0 n/ _        z=z+s*y;
      }
    return(z);
  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
  int n,m;
  double x[],y[],a[],dt[];
3 z# ^8 p3 _: C$ @$ v  { int i,j,k;
& X* w; g7 V# F0 S5 e! e    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
" \; Z' b7 I/ Y; O3 f- S    for (i=0; i&lt;=m-1; i++) a=0.0;
- \' @' W& d! k8 @  U# A    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);
6 z# A% d. m! O" a# R    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;
8 x& h3 F4 k7 p' N* e: v6 |    a[0]=c*b[0];
    if (m&gt;1)
4 C" T, N, D% a- F      { t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
6 f& k! n0 H8 f5 U6 k3 J! I+ l        for (i=0; i&lt;=n-1; i++)
          { q=x-z-p; d2=d2+q*q;
/ @# p) v& T# j3 {1 [3 H' P            c=c+y*q;
            g=g+(x-z)*q*q;
          }
        c=c/d2; p=g/d2; q=d2/d1;
5 N- g" p6 d# U. ]" V        d1=d2;
& v( W$ F2 R. y4 `$ r        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
    for (j=2; j&lt;=m-1; j++)
      { s[j]=t[j-1];
' a8 {9 A5 {1 p0 N        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
% q/ ]% d+ Z* N; ?          for (k=j-2; k&gt;=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
' E' |- H6 X/ @) a        s[0]=-p*t[0]-q*b[0];
5 K' n+ Z" h7 x2 v0 H        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
: B9 }; b! B! J& d! v( L8 A& r            for (k=j-1; k&gt;=0; k--)
              q=q*(x-z)+s[k];
3 I0 D' }4 z. }% c5 ?# e            d2=d2+q*q; c=c+y*q;
/ H' |7 l1 i. y3 f. w            g=g+(x-z)*q*q;
7 }" s' _; V! P8 c          }
        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];
            b[k]=t[k]; t[k]=s[k];
          }
      }
* g/ L: O+ o) w  o6 ]. [; m/ I' ?    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
- j" C' N8 F6 J6 y& R) N/ w; }    for (i=0; i&lt;=n-1; i++)
      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
" I8 k: _9 ?+ ^, e* P+ V5 X7 ~' a          q=a[k]+q*(x-z);
        p=q-y;
2 ^$ Y$ W; a# I( a1 D) G        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
$ S2 J0 Y1 H% z/ A2 x2 l        dt[0]=dt[0]+p*p;
2 G. g: w( O2 f4 M! Y$ L% K/ X        dt[1]=dt[1]+fabs(p);
      }
    return;
6 T& I( J" n7 r6 R1 w* H  }</P><>龙贝格积分法</P><>#include "math.h"
: q# {6 _( y  r4 q) a+ E8 `  q0 L  double fromb(a,b,eps)
2 r# N8 o  E5 u6 e  double a,b,eps;
* J' |  C" B/ Z& L7 q  u) w  { extern double frombf();
    int m,n,i,k;
0 j$ K- Z3 ~. e3 |: z: S6 y% y. X    double y[10],h,ep,p,x,s,q;
    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
: _# N6 B. F! g. H    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
      { p=0.0;
        for (i=0;i&lt;=n-1;i++)
0 Q, r8 Z% x0 b, B8 j- H          { x=a+(i+0.5)*h;
            p=p+frombf(x);
  S' X9 _" k# z$ g% u          }
: x4 E# m; q! q& t+ D        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k&lt;=m;k++)
          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
' r3 h6 I( B1 ~- K            y[k-1]=p; p=q;
- A$ j6 B  X' l* b# y+ H          }
4 C  k# S) S3 F8 a2 W, j2 C4 {1 S        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
/ `/ M+ |# C# I/ M" @      }
    return(q);
  }</P><>呵呵 希望对你有用！！</P>

西西

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

陈了

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