数学建模社区-数学中国

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

---

作者: 西西    时间: 2004-5-10 08:46
标题: 求助！！ 数值分析程序（c语言）
<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
: f; b0 F, R( A3 R<>1。牛顿迭代法</P>
# U8 V  i; [* T. `<>1。全主消元法</P>
" z2 r* P& t3 }& l5 E/ P<>3。改进平方根法</P>
8 o5 Q6 j. m; q# H) j9 e7 l8 S<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
5 Q: _$ H7 Y; J5 C/ I- |  G% e; Y0 Z<>6。龙贝格积分法</P>
, a3 F* y& n$ J- R<>拜托了。。。。</P>

---

作者: loveyuan    时间: 2004-5-11 08:32
<>1。牛顿迭代法</P><>
  #include "stdio.h"
  #include "math.h"
7 S3 i3 y! n# m# U  int dnewt(x,eps,js)
3 [( s# `, y8 y! m  int js;
5 s3 g( c5 b* Y) ^  double *x,eps;
4 E1 {, F- u+ v1 |5 s' m+ n  { extern void dnewtf();
/ G8 [3 M! b- }3 ?/ e+ b    int k,l;
    double y[2],d,p,x0,x1;
# S8 ]4 j2 N' y$ _! C5 m8 ]. N    l=js; k=1; x0=*x;
    dnewtf(x0,y);
    d=eps+1.0;
    while ((d&gt;=eps)&amp;&amp;(l!=0))
      { if (fabs(y[1])+1.0==1.0)
& H( F' d6 Q2 k- k          { printf("err\n"); return(-1);}
        x1=x0-y[0]/y[1];
        dnewtf(x1,y);
7 f8 \' y4 s& D9 A        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
3 o) j9 l1 K) M* ~+ i+ L4 Q$ ~. Y      }
    *x=x1;
    k=js-l;
( ], p1 n; L. Z# M+ O+ K: K, l    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
( }6 \7 u7 ]+ ~5 t; i& j2 E  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
' o# t) V  w# M; E( E# _    double p,q,s,d;
3 F8 Y- E" B% M' w8 o3 M) s    js=malloc(n*sizeof(int));
, n' t# t3 u3 y3 V4 r* z    for (k=0;k&lt;=n-2;k++)
; J3 R7 F7 }) o4 i% M# w- p' Z      { d=0.0;
& T& z- _- D' a" p! i, D        for (i=k;i&lt;=n-1;i++)
- k, G; V5 j& S: e        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;}
8 q( O  p. Y0 l* m          }
9 w' a4 Q) z' N2 ~* y  z% L4 D! y6 h+ b        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
            return(0);
& Y$ |( d6 F. [0 M+ P+ |          }
0 d+ \+ W+ T& F3 F2 `+ d        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
+ ~/ H  F4 N3 V- N$ j& |1 D# m  i              { u=k*n+j; v=is*n+j;
% g' z  c# f. E                p=ar; ar=ar[v]; ar[v]=p;
% o3 ~& X1 h* ]7 w% C2 N7 F' ~2 }& ?                p=ai; ai=ai[v]; ai[v]=p;
              }
% ?) ^4 m2 c* I  F! l            p=br[k]; br[k]=br[is]; br[is]=p;
6 s2 B! u/ r; i' M6 W0 F4 n) q            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
# ~3 @% w' {- ?8 ^& t$ d0 A          }
, _" T  W% q" f' k' l! G+ c        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
8 }+ v  c, e& ?5 t8 m9 h( y            { 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;
            }
. D, G) ^' _* a$ g( Y/ N        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
' P" S) r" Q7 {          { u=k*n+j;
/ u. m9 K) L: g% }# \1 X# @- V            p=ar*ar[v]; q=-ai*ai[v];
            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
" f1 {$ ^6 m& D5 e" {8 x! B        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i&lt;=n-1;i++)
- r$ [+ N, c8 S" E          { u=i*n+k;
            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];
" h9 K0 w7 X9 Q/ @, v4 t8 p                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
3 O. H- T( b; `2 T            p=ar*br[k]; q=ai*bi[k];
$ [9 S+ {- e% W! `; t7 K6 e* _            s=(ar+ai)*(br[k]+bi[k]);
            br=br-p+q; bi=bi-s+p+q;
          }
( \# R) e1 r. _      }
    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
: D& n3 X& x( N0 [% @( V        return(0);
9 Y$ R2 L0 i6 R; n      }
) f, p, |$ l' j- `    p=ar*br[n-1]; q=-ai*bi[n-1];
7 P9 [7 b2 I2 \9 L0 v    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--)
7 d! f9 H! ]( [+ |8 _$ F; Y    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
$ [! Z' c& u+ G6 A  u9 d4 U        p=ar*br[j]; q=ai*bi[j];
3 P  {8 z. T  Q        s=(ar+ai)*(br[j]+bi[j]);
) ~% y; O) w& w+ I- A        br=br-p+q;
        bi=bi-s+p+q;
0 e3 e  t' ]4 a' L$ O      }
/ I2 h. h3 \1 j8 b) R2 h    js[n-1]=n-1;
) M2 ]$ f0 q3 [) i# e% @7 c    for (k=n-1;k&gt;=0;k--)
' i4 F0 l6 [4 W; }# d      if (js[k]!=k)
8 l" J) v3 L0 B# ?5 |) M9 a        { 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);
" A* U; R1 ]) S$ R& G' U9 U    return(1);
  }</P><>平方根法</P><>#include "math.h"
  #include "stdio.h"
, X' t$ [( }+ R0 g- L+ r8 a1 w/ }  int achol(a,n,m,d)
  int n,m;
$ E6 R4 f. p8 e  double a[],d[];
  { int i,j,k,u,v;
) {* {7 Q) Z1 v2 P7 V    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
# a% A- a- ~- R$ k+ ^. F      { printf("fail\n"); return(-2);}
- J; l. W5 A& B2 o( o- o$ \) b    a[0]=sqrt(a[0]);
- ]- N$ L8 U+ E$ R    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
+ j8 ~5 c5 _: p8 R5 T- `      { u=i*n+i;
- E" [% s6 D& V3 F5 f5 |) a        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
; [1 R# R# H! S& D          }
        if ((a+1.0==1.0)||(a&lt;0.0))
) h: ]: o2 m# r4 _6 E          { printf("fail\n"); return(-2);}
        a=sqrt(a);
        if (i!=(n-1))
          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
) Q( _( H+ s. F# k3 ~. j# S' J                for (k=1; k&lt;=i; k++)
1 |/ e' W5 b# i* h9 Y$ y: _. y6 x. A                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
* o) ?4 n0 M' d1 }8 x' i& r5 w                a[v]=a[v]/a;
              }
          }
      }
    for (j=0; j&lt;=m-1; j++)
      { d[j]=d[j]/a[0];
8 B( b' e5 e. x  P: ~/ i5 y        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
1 u( r; X6 y. U! D            for (k=1; k&lt;=i; k++)
3 X  O8 r, d/ A# Y' A              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
            d[v]=d[v]/a;
6 N3 W3 S/ [6 |6 B          }
# t0 Q# C9 `" Q      }
6 B( |5 _5 g  f* h7 c    for (j=0; j&lt;=m-1; j++)
# f1 y3 Y5 o8 c      { u=(n-1)*m+j;
# X6 Y) `! t, d2 ?3 l' N: r0 U        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++)
% i: t8 A  Q/ _" |6 @              { v=(k-1)*n+i;
, N* R. k7 `+ ]( Q8 J. c! w                d=d-a[v]*d[i*m+j];
+ C9 p1 |" Z" n9 `/ e' ~" y              }
4 {3 b8 l' `) ~' v8 @7 _            v=(k-1)*n+k-1;
            d=d/a[v];
" a8 O8 C$ c7 S. F+ |5 Y  K          }
2 }- G$ c6 P6 Y- `0 [  I' ]      }
    return(2);
  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
  { int i,j,k,m;
, |/ v9 \7 @) R! |! K! g4 d    double z,s,xi,xj;
    float p,q;
9 L8 s- c  A1 ~: i    z=0.0;
" o8 g; o& p: K* p% v5 R    if (n&lt;1) return(z);
1 z, J* X" b/ ]9 I' r& s    if (n==1) { z=y[0]; return(z);}
    if (n==2)
# M7 O" M/ A' n! a/ B+ W" y      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
: e+ ?' h/ Z5 Y/ g9 I9 X' e& |* y- m        return(z);
      }
    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
, S# A! ~: ?3 w        if (p&gt;q) i=i+1;
      }
% o) O2 H0 j, o1 }: u. O7 p    else i=0;
$ B$ {8 v; N* E( {    k=i-4;
% E* ~* x& v- `7 ]$ _    if (k&lt;0) k=0;
9 g; w  t6 r3 P# m    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++)
# y; p5 ^) R4 [7 w          if (j!=i)
            { xj=x0+j*h;
$ f2 x: k9 b/ }, P& V0 d0 @% p              s=s*(t-xj)/(xi-xj);
/ J6 @3 [8 C  ?, V& Z) h3 m            }
" h, D" q9 E$ V7 }9 k        z=z+s*y;
      }
5 s6 _+ w6 `6 Y7 J7 Z    return(z);
  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
1 z% W, `& B% D* V6 R  k  void hpir1(x,y,n,a,m,dt)
; Q8 v6 C( x$ K  ^6 z0 L6 R3 @! X  int n,m;
( B4 @+ C* n0 V" P$ v  double x[],y[],a[],dt[];
# X2 z* {$ e7 {! a  { int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
6 r0 }# B. I9 ~% G3 Y/ K1 G# \5 m) w/ g8 ?    for (i=0; i&lt;=m-1; i++) a=0.0;
    if (m&gt;n) m=n;
+ e$ \0 e) `# _: i$ j1 n    if (m&gt;20) m=20;
    z=0.0;
: L$ U0 R" B: r+ z    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
) K+ }: W* w1 s2 s) A    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
* ]; G3 J2 b1 B- t) y( d: w    for (i=0; i&lt;=n-1; i++)
9 u$ ?- k$ E, d" Y1 D4 x7 |      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
2 ^' W- q1 S( f" S    a[0]=c*b[0];
    if (m&gt;1)
% w' h) j! j7 _6 F      { t[1]=1.0; t[0]=-p;
5 m- b; m2 d1 S$ b: @6 b        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
. m) E: I! P& I. w* \          { q=x-z-p; d2=d2+q*q;
: S* i. j* B& D1 f            c=c+y*q;
            g=g+(x-z)*q*q;
2 F6 \& f( r+ _0 r& `6 N0 G          }
) }2 A, ~/ c: F8 f& K        c=c/d2; p=g/d2; q=d2/d1;
5 x- S) B: {, p' H" T% Z& y        d1=d2;
4 v1 i3 [/ o+ X8 o" _& u1 f        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
/ x  _6 U' B) J& r7 ]      }
' Q5 h$ U" \/ g- n* h, `: T    for (j=2; j&lt;=m-1; j++)
      { s[j]=t[j-1];
        s[j-1]=-p*t[j-1]+t[j-2];
( H7 t# j) V+ y% q$ P        if (j&gt;=3)
          for (k=j-2; k&gt;=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
        s[0]=-p*t[0]-q*b[0];
1 J! X- q/ D5 b& \9 w        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
0 Q* ^) S. g' T4 f1 @: }          { q=s[j];
            for (k=j-1; k&gt;=0; k--)
( z& o1 n4 P4 ~! |              q=q*(x-z)+s[k];
& M+ K9 w) t0 L) ^, f3 k            d2=d2+q*q; c=c+y*q;
, j) }: Q0 l& a            g=g+(x-z)*q*q;
, a5 Q" }; |8 M2 I' j, l          }
        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--)
. S. s4 Z: z+ `& ^          { a[k]=c*s[k]+a[k];
7 e" `) s. C* ~5 K% U            b[k]=t[k]; t[k]=s[k];
. p" N* J4 Z8 r( ?6 g* e          }
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
    for (i=0; i&lt;=n-1; i++)
) Q+ [1 ~& e9 }7 s! z( G" W" k      { q=a[m-1];
  }6 O9 q5 Y9 J        for (k=m-2; k&gt;=0; k--)
          q=a[k]+q*(x-z);
        p=q-y;
* U' z( K: k( r  e, x/ m& T        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
# }1 t6 `3 t8 }* h) ?" N        dt[0]=dt[0]+p*p;
- P' W) M" p- E' x1 s        dt[1]=dt[1]+fabs(p);
      }
    return;
6 Z- f- B* m8 A3 O0 X  }</P><>龙贝格积分法</P><>#include "math.h"
7 z3 T3 ?7 _1 z, S7 W  double fromb(a,b,eps)
  double a,b,eps;
; m5 Z# Z9 [! D: D; A- n  { extern double frombf();
    int m,n,i,k;
+ D% ^8 m7 D; t* {0 c: J    double y[10],h,ep,p,x,s,q;
' x5 v* A0 C8 h4 X    h=b-a;
+ A/ ]' S+ F+ {6 R    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
      { p=0.0;
        for (i=0;i&lt;=n-1;i++)
# j# T4 t) |0 t          { x=a+(i+0.5)*h;
0 j3 c9 v+ N% k3 u7 i9 k            p=p+frombf(x);
- {! n# C1 R' U/ @; u8 O* D* C. _          }
        p=(y[0]+h*p)/2.0;
5 ^0 ]& @- ~3 n  \- b- u) C' y        s=1.0;
        for (k=1;k&lt;=m;k++)
          { s=4.0*s;
! ^  V$ R4 `! S            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
# T5 b+ i. R1 N- l      }
    return(q);
% y4 N$ {9 ?; _4 T  }</P><>呵呵 希望对你有用！！</P>

---

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

---

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

---

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

Powered by Discuz! X2.5