数学建模社区-数学中国

标题: 求助！！ 数值分析程序（c语言） [打印本页]

---

作者: 西西    时间: 2004-5-10 08:46
标题: 求助！！ 数值分析程序（c语言）
<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
5 M4 M, L. u, Y<>1。牛顿迭代法</P>
. v2 l0 B9 N: g4 _! g: t: a, A<>1。全主消元法</P>
( U- w0 ^" `: q<>3。改进平方根法</P>
) f. N5 w* g# @) C  Q* J' Q<>4。牛顿向前，向后插值法</P>
4 m+ v; y* x: L% H$ D. }5 U' O<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
<>拜托了。。。。</P>

---

作者: loveyuan    时间: 2004-5-11 08:32
<>1。牛顿迭代法</P><>
' J6 R; Q. l$ j  j0 t  #include "stdio.h"
  #include "math.h"
( ~7 |$ H3 h- Q" r  int dnewt(x,eps,js)
  int js;
  double *x,eps;
0 {5 G& {. M* `4 _  { extern void dnewtf();
, B8 c( D) w8 C1 {# Q    int k,l;
    double y[2],d,p,x0,x1;
! V8 L0 x/ K, k. u9 k7 {0 T/ ^( A8 W% ~    l=js; k=1; x0=*x;
    dnewtf(x0,y);
    d=eps+1.0;
/ k2 A* X1 @4 K* D    while ((d&gt;=eps)&amp;&amp;(l!=0))
      { if (fabs(y[1])+1.0==1.0)
          { printf("err\n"); return(-1);}
( I7 p  B7 w* A$ O, K. a" x# v        x1=x0-y[0]/y[1];
        dnewtf(x1,y);
        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
      }
# ~* S, g1 f' {6 M( b    *x=x1;
- e2 z! P4 @: m0 ]    k=js-l;
1 ]0 g# W9 ^3 e+ D+ ]    return(k);
5 _# e5 }- g0 v8 p( W% ?' R% p  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
, `4 L4 x  Y! T. p- l  int acgas(ar,ai,n,br,bi)
, P! m$ v) B' r6 G- N1 K& x; F  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
* {6 r- ~- i" v+ D    double p,q,s,d;
    js=malloc(n*sizeof(int));
( ^, ]# _' J# S4 S1 G    for (k=0;k&lt;=n-2;k++)
! T9 o# y" M( Z- H' J      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
" |% V( \& @$ B  L" f- P* O4 w        for (j=k;j&lt;=n-1;j++)
5 x6 Z  ^+ _! I9 Q7 F2 @. K: p- w& P" S          { u=i*n+j;
- z& r! k- L4 u2 i( ~            p=ar*ar+ai*ai;
/ u* x: a0 A& @. F% X6 Q6 Z            if (p&gt;d) {d=p;js[k]=j;is=i;}
          }
: D8 u, R+ l+ m4 O        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
  |4 \3 ~+ w0 d$ ?  ]            return(0);
5 m7 U7 g6 E' b# D6 Y0 h; B) U7 B$ w          }
        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
6 w$ V4 h1 Y) u8 Q$ Z; u& d2 F                p=ar; ar=ar[v]; ar[v]=p;
* ~0 I; r& ~; V8 K                p=ai; ai=ai[v]; ai[v]=p;
  g8 c0 g) N2 |! G7 Z3 l1 v              }
5 l# C4 k" Y! m  I9 B* A7 E            p=br[k]; br[k]=br[is]; br[is]=p;
            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
          }
        if (js[k]!=k)
          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;
8 m! k, @0 i5 f8 {) P2 u  X$ _            }
        v=k*n+k;
9 W2 t4 j$ V4 J9 O' Z  l5 ?% w        for (j=k+1;j&lt;=n-1;j++)
- P! h0 @* P! h2 u4 V. \          { u=k*n+j;
% J) k: ]& b+ T: s            p=ar*ar[v]; q=-ai*ai[v];
            s=(ar[v]-ai[v])*(ar+ai);
; k/ c/ b6 W5 H/ F            ar=(p-q)/d; ai=(s-p-q)/d;
5 w% Q4 E  q1 ?          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
! u( U1 R& S6 u; f, G, n8 n# r        s=(ar[v]-ai[v])*(br[k]+bi[k]);
& g. R0 e: l3 [: T5 ?/ N2 R, ]5 F% G        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
& y$ ]& e, m! r. F' j& d        for (i=k+1;i&lt;=n-1;i++)
2 a: y- [/ H0 H6 F: ~& A, ~          { u=i*n+k;
& e4 F- Z6 k' @( F* ~6 \            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]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
! J2 T2 |' ]$ E+ W- m& ^              }
5 H9 a# ^" i: [3 G            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
+ L1 G- l0 f8 ~/ v  v0 m' q            br=br-p+q; bi=bi-s+p+q;
          }
      }
    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
; c$ {% u$ d0 n3 ?      { free(js); printf("err**fail\n");
        return(0);
  Y$ g! v% g/ Z$ e" K* G3 Y* C$ A      }
    p=ar*br[n-1]; q=-ai*bi[n-1];
  d% s& L" D' T1 d; |( V    s=(ar-ai)*(br[n-1]+bi[n-1]);
    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
" ^) x# |" A" i- y9 C    for (i=n-2;i&gt;=0;i--)
& V$ w1 h$ c) t# d- m" Q    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
! I( s, V% y# G+ u/ P        p=ar*br[j]; q=ai*bi[j];
: Z: C1 |2 j: G. ^+ u7 U        s=(ar+ai)*(br[j]+bi[j]);
. R# o. h. R0 H. [        br=br-p+q;
        bi=bi-s+p+q;
      }
    js[n-1]=n-1;
% a9 j4 T2 G  Y; f0 }( ]# i5 m/ L    for (k=n-1;k&gt;=0;k--)
5 `: _3 |4 t$ x( O% L4 H      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;
        }
    free(js);
    return(1);
, C- E4 |( c$ k1 r4 n" T* C8 ]  }</P><>平方根法</P><>#include "math.h"
/ p3 I& b0 H/ g8 \) e  |' H( n  #include "stdio.h"
  int achol(a,n,m,d)
2 g" v+ S) z+ [6 e/ ^) ^  int n,m;
  double a[],d[];
  { int i,j,k,u,v;
    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
: n- `$ L4 j7 i$ x) l      { printf("fail\n"); return(-2);}
" `  [0 p4 n! C+ p  P; |- P3 V* m    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++)
4 e  H5 C. O  Q6 H      { u=i*n+i;
        for (j=1; j&lt;=i; j++)
3 {" D( f3 V; Z" ^1 f$ L! E' y          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
          }
$ u8 i- k/ T5 G' r0 Y        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
* J  U& ?, K. G4 Y) j        a=sqrt(a);
- q& j+ F8 {3 b6 O6 F        if (i!=(n-1))
          { for (j=i+1; j&lt;=n-1; j++)
7 Q- i& u( y  J/ o" P9 ~. f! X              { v=i*n+j;
# g& |# M; |! l# J# ~* ~3 {# m( N0 C* A                for (k=1; k&lt;=i; k++)
# m) N% a& D2 n1 x0 g                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
2 N0 p* B4 N; L/ @) C3 g0 c                a[v]=a[v]/a;
              }
          }
9 j' Q' X3 J- B# v4 z6 ^      }
  r' _. I. s6 z' T7 X3 C    for (j=0; j&lt;=m-1; j++)
2 ]: |; z3 I8 g4 S      { d[j]=d[j]/a[0];
        for (i=1; i&lt;=n-1; i++)
4 ?# ]! |, o% b0 Q          { u=i*n+i; v=i*m+j;
            for (k=1; k&lt;=i; k++)
$ K9 }+ I: H" |: X$ ~              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a;
2 R3 R- _7 T$ z1 b          }
4 Q! `2 P- R$ |% K, X4 Q      }
    for (j=0; j&lt;=m-1; j++)
      { u=(n-1)*m+j;
# j9 Z0 p: B0 o" |7 b. g% C: Y* p        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++)
  a/ S0 Q5 f9 M0 X( K, Q3 k* Q2 U              { v=(k-1)*n+i;
5 E  s" \& |% \- p                d=d-a[v]*d[i*m+j];
& @- Q) ]' t% n              }
            v=(k-1)*n+k-1;
9 _( L  Q0 f4 ]  G( A& r            d=d/a[v];
: n+ v9 j/ o. e+ ]. Q% R          }
      }
    return(2);
4 G4 g/ M0 u+ z' |' J3 Q- ?  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
; Z) Q3 `5 Q: A6 Y8 p/ E  { int i,j,k,m;
    double z,s,xi,xj;
& D' g6 Y* Q# `, x- g8 [    float p,q;
    z=0.0;
    if (n&lt;1) return(z);
$ m& V' B5 ~/ w    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
( ^( Q' s5 ~/ @/ v: T* O        return(z);
6 S3 e' Y7 o: Z' }6 p/ E      }
    if (t&gt;x0)
+ W9 x5 _. l( J" q$ a- W2 @- p      { p=(t-x0)/h; i=(int)p; q=(float)i;
8 U  b/ I, G* d1 q- r        if (p&gt;q) i=i+1;
      }
9 g8 w8 W, Z9 L0 v    else i=0;
: b9 Z( H5 q6 F, z6 _7 l    k=i-4;
    if (k&lt;0) k=0;
2 S+ `+ r, b& q6 ?% J$ i0 f    m=i+3;
9 `" M  {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++)
8 M% [% _8 S/ ]- @# P! ^. k          if (j!=i)
            { xj=x0+j*h;
3 K0 a. C5 w9 K1 _: R4 }! r              s=s*(t-xj)/(xi-xj);
            }
8 i+ X& ^( \$ A+ T9 B        z=z+s*y;
: U  p$ D; W! g( t5 l1 v7 X9 J) z' _      }
7 {" a" X! M* w- B" q1 J    return(z);
+ L' g+ q1 X1 X0 d3 ?  }
6 a3 ^! X# s! i. q% ]1 K8 b2 y向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
: P3 ^% F& d; O! N8 H8 f5 d  void hpir1(x,y,n,a,m,dt)
+ e: p: L) k4 e4 N7 Y  int n,m;
2 g3 ?. d( Q+ _! l0 p* n* \  double x[],y[],a[],dt[];
! _7 Q& q: w' E' f; V8 J' I9 s' i. m1 e  { int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
; s+ _# I0 j* }- |    for (i=0; i&lt;=m-1; i++) a=0.0;
& @2 E3 F1 c* F6 _' H- k+ j0 c    if (m&gt;n) m=n;
7 o1 B0 E1 d7 p    if (m&gt;20) m=20;
    z=0.0;
& j- a4 A; X. ?( P/ \8 ^    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
! f2 t: l9 o4 I- k+ H: @* \    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
- r; I4 |/ T7 |    for (i=0; i&lt;=n-1; i++)
5 Q+ p3 K* b7 `. ]" @4 v      { p=p+(x-z); c=c+y;}
; j7 _* q. K2 Y  v' `    c=c/d1; p=p/d1;
6 u: @, T  E' v9 V    a[0]=c*b[0];
. W  P  F, {5 {! d6 E    if (m&gt;1)
      { t[1]=1.0; t[0]=-p;
/ B4 M/ }1 p( b( N5 e& a3 q! i. ~        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
9 m! C7 f1 k# G6 `- i          { q=x-z-p; d2=d2+q*q;
            c=c+y*q;
" L  G, A  s/ g7 P; E            g=g+(x-z)*q*q;
# e0 G. d' f% L) U3 H2 Y/ V2 F          }
8 [/ J) m3 r* j& k4 H" i9 U' U        c=c/d2; p=g/d2; q=d2/d1;
        d1=d2;
  T% P+ ~/ }+ I/ u        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
    for (j=2; j&lt;=m-1; j++)
8 [( x; Z$ S  V& G- [9 U0 O) G0 i, q      { s[j]=t[j-1];
- B* r, w& ^2 A/ e- e        s[j-1]=-p*t[j-1]+t[j-2];
" S7 X1 S) |3 h9 b* o        if (j&gt;=3)
          for (k=j-2; k&gt;=1; k--)
) Q5 l* ]$ V% E2 v            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++)
; @7 w+ [2 C1 d0 r/ z          { q=s[j];
& Q+ _( T0 h% z* k- ^$ K! c* r            for (k=j-1; k&gt;=0; k--)
5 y) N; _, q- i8 y: [. C% ^              q=q*(x-z)+s[k];
( W1 \$ L% x. V3 B            d2=d2+q*q; c=c+y*q;
1 x: m$ J- ^; V9 V8 M- z5 R            g=g+(x-z)*q*q;
          }
        c=c/d2; p=g/d2; q=d2/d1;
3 [4 ~/ h1 ?8 E- c$ i) s$ V        d1=d2;
$ ^: w7 v2 ?+ k3 u& G8 ^: D+ V* H! v/ w        a[j]=c*s[j]; t[j]=s[j];
        for (k=j-1; k&gt;=0; k--)
+ K! m$ c' g; o& O0 j          { a[k]=c*s[k]+a[k];
            b[k]=t[k]; t[k]=s[k];
          }
      }
    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--)
/ |+ z$ v/ ~0 t6 W0 f7 N+ |          q=a[k]+q*(x-z);
        p=q-y;
6 f2 l+ \) N- e/ L" S- q# I        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
6 n$ U4 X( {/ [; l4 y% E        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
. Z( m7 c& t% p3 n9 {      }
    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
9 {7 q& c* P* ^" i7 Q  double fromb(a,b,eps)
, e+ r4 }2 n7 G, g7 }. J; X. D  double a,b,eps;
  { extern double frombf();
2 X  R5 P9 X- l( [  j9 T/ {" F    int m,n,i,k;
    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;
* q, s: g. y; s$ e/ O    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
; x( @) a: g  ~) s0 H      { p=0.0;
+ c: ]- N+ J! V* c        for (i=0;i&lt;=n-1;i++)
; M5 f" H5 r# w, u- R1 R          { x=a+(i+0.5)*h;
            p=p+frombf(x);
          }
        p=(y[0]+h*p)/2.0;
2 I' x5 R& Y7 t        s=1.0;
1 @3 i/ {: l6 k/ c% W        for (k=1;k&lt;=m;k++)
          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
2 t4 a4 V0 [( y5 _' r) c) ?            y[k-1]=p; p=q;
1 b- B6 F7 O/ [7 S) X9 s; k          }
0 B0 g. \) y4 o" t0 a        ep=fabs(q-y[m-1]);
8 Q4 O3 b) W$ z% E9 v6 G& X        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
: x% ?* d+ e. K8 `" \; L  _% C  s. @      }
    return(q);
  }</P><>呵呵 希望对你有用！！</P>

---

作者: 西西    时间: 2004-5-14 10:54
标题: ？
<>谢谢你哦，呵呵</P>

---

作者: 陈了    时间: 2010-3-22 11:13
哇哇哇哇哇哇。。。好复杂呀！看不懂

---

欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)

Powered by Discuz! X2.5