数学建模社区-数学中国

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

---

作者: 西西    时间: 2004-5-10 08:46
标题: 求助！！ 数值分析程序（c语言）
<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
) Y5 h- o  Q7 Q9 d* Q6 S" i6 g: _; N<>1。牛顿迭代法</P>
: Q# c* d6 n+ P! s6 j4 I$ a<>1。全主消元法</P>
<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
/ V! v! p9 }- K. N  e/ I+ E/ ?0 w<>拜托了。。。。</P>

---

作者: loveyuan    时间: 2004-5-11 08:32
<>1。牛顿迭代法</P><>
1 ]6 S1 j9 F9 s/ N  #include "stdio.h"
  #include "math.h"
  int dnewt(x,eps,js)
  int js;
, b4 ?1 X& ]' J) G' q8 {& ~# u  double *x,eps;
  { extern void dnewtf();
" j9 F  s  o- p! a3 s( r  i& T    int k,l;
9 l( N3 g" R, f    double y[2],d,p,x0,x1;
( r- d6 ?8 I( ^, B    l=js; k=1; x0=*x;
    dnewtf(x0,y);
  }0 O6 i; s! L( Q: ^  E    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];
1 h1 B+ D6 s; q6 S0 [$ i* O5 C        dnewtf(x1,y);
. F* Q' `8 J+ h4 s5 `# A        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
7 ^7 E- T- z  ]+ g5 O% T3 j        x0=x1; l=l-1;
" H% i" G3 ^- P+ M' J      }
* U1 v4 u5 m3 r    *x=x1;
    k=js-l;
& M' ^* c: \# ^/ c5 v    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
  int n;
  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
    double p,q,s,d;
    js=malloc(n*sizeof(int));
1 v" h& V  y7 K/ x    for (k=0;k&lt;=n-2;k++)
9 H1 |$ B" [, f$ U      { d=0.0;
- I" |1 V9 @" U7 e        for (i=k;i&lt;=n-1;i++)
        for (j=k;j&lt;=n-1;j++)
          { u=i*n+j;
            p=ar*ar+ai*ai;
' r% V, f2 }$ H/ n) x            if (p&gt;d) {d=p;js[k]=j;is=i;}
& d0 L. K. ?  r! f          }
        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
            return(0);
" _$ p6 q$ c9 P( j, O          }
        if (is!=k)
          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
- @# ?, h6 y0 E( O, V                p=ar; ar=ar[v]; ar[v]=p;
                p=ai; ai=ai[v]; ai[v]=p;
) P2 u/ \+ U4 C: H  X              }
: R6 g. f4 O/ }7 m" }/ w            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++)
; D6 J: }7 L+ D7 o$ ?! Y  y; A            { u=i*n+k; v=i*n+js[k];
              p=ar; ar=ar[v]; ar[v]=p;
9 v; _" B+ C+ X8 b              p=ai; ai=ai[v]; ai[v]=p;
- L1 t2 u# O4 E. k+ E' D/ r            }
5 ~* v5 R: j+ Z: B( X/ L        v=k*n+k;
  p6 ~1 w8 Z4 k        for (j=k+1;j&lt;=n-1;j++)
0 X8 k7 d3 Y- a% g3 j, U          { u=k*n+j;
            p=ar*ar[v]; q=-ai*ai[v];
  ?6 D: b+ X: s0 Q7 j# M- V            s=(ar[v]-ai[v])*(ar+ai);
            ar=(p-q)/d; ai=(s-p-q)/d;
          }
% T7 t' [; {! M1 B7 X        p=br[k]*ar[v]; q=-bi[k]*ai[v];
        s=(ar[v]-ai[v])*(br[k]+bi[k]);
# k& q! n1 f: Q        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i&lt;=n-1;i++)
$ _5 b4 J; h( d" G( n          { u=i*n+k;
            for (j=k+1;j&lt;=n-1;j++)
; i8 i" |4 t  N/ _  A              { v=k*n+j; l=i*n+j;
                p=ar*ar[v]; q=ai*ai[v];
  t9 N1 E5 E4 l! ~4 M9 k  R% g                s=(ar+ai)*(ar[v]+ai[v]);
) t  [4 ^: z/ w3 h2 a- _1 w                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
$ M$ O/ b1 Y0 C2 K* d            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
) \( |* D; A( [6 I) D            br=br-p+q; bi=bi-s+p+q;
, J: x" O4 S9 s) ^2 m. J          }
+ V  H: F9 i3 U$ D      }
! W; I* f2 c' ]    u=(n-1)*n+n-1;
    d=ar*ar+ai*ai;
    if (d+1.0==1.0)
9 g1 {$ B) T" V      { free(js); printf("err**fail\n");
        return(0);
      }
# ]8 A2 `9 U; }# [- K    p=ar*br[n-1]; q=-ai*bi[n-1];
; m. x6 b$ `0 V7 S+ ~" q3 d' d  Y8 S    s=(ar-ai)*(br[n-1]+bi[n-1]);
+ P# J2 G1 O, l/ e& T  V    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
8 P) Z0 B3 F% C8 _8 ]) i) {% b    for (i=n-2;i&gt;=0;i--)
    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
- U+ |  V- P& _" K1 e        p=ar*br[j]; q=ai*bi[j];
' x. N3 U0 N( L7 X        s=(ar+ai)*(br[j]+bi[j]);
        br=br-p+q;
        bi=bi-s+p+q;
      }
    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;
          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
7 F& Q/ l7 h7 `  W        }
    free(js);
    return(1);
4 u" u8 s. G% H) ~, t& g( s$ L  }</P><>平方根法</P><>#include "math.h"
: c8 l8 a6 \7 o  #include "stdio.h"
+ M' M2 t" m6 s; l& O  int achol(a,n,m,d)
- Y* o" o7 ^7 W! X) i  int n,m;
  double a[],d[];
! c# [% N' g4 V( W  { int i,j,k,u,v;
( S& T, @* |. b; w; U    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
) H7 _( |3 f/ P0 b1 N' Y. Q+ z. D! n      { printf("fail\n"); return(-2);}
. m+ e; u3 M8 j" D    a[0]=sqrt(a[0]);
$ v4 t' M3 c% G/ x: E' i, l" K5 G/ p    for (j=1; j&lt;=n-1; j++) a[j]=a[j]/a[0];
    for (i=1; i&lt;=n-1; i++)
      { u=i*n+i;
) i! u* Q# s6 @# j! t        for (j=1; j&lt;=i; j++)
) G: V( s1 \  J/ \          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
8 {  ~' i% Y/ N( n7 X9 P% Q  |/ z          }
        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
$ d4 X2 I( g: _7 N! l+ v        a=sqrt(a);
        if (i!=(n-1))
4 B0 F* X" |, C0 u+ D& L! {% N          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
                for (k=1; k&lt;=i; k++)
2 q& R6 ~4 x5 W                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
                a[v]=a[v]/a;
8 X& j3 j3 M& [( h7 z4 B6 S              }
( A5 O8 j4 |$ C9 b/ {          }
* ^  _) E. m1 @0 [      }
    for (j=0; j&lt;=m-1; j++)
      { d[j]=d[j]/a[0];
6 A( u' p1 [4 j; y8 {4 z2 g/ l* M        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
; W; l) p( }% @, U9 P            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;
          }
      }
; S! s. W6 l0 x    for (j=0; j&lt;=m-1; j++)
" I0 @. _( t, r      { u=(n-1)*m+j;
: y# s6 d9 a: U7 Z& W$ G# j        d=d/a[n*n-1];
        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
# @8 r& C5 i; b- V$ Z            for (i=k; i&lt;=n-1; i++)
              { v=(k-1)*n+i;
: }) Y( \9 e0 r" n+ V+ R: q                d=d-a[v]*d[i*m+j];
              }
2 B1 K$ p3 I  p# x4 A+ V6 c5 l            v=(k-1)*n+k-1;
$ v: K7 Y% l' ^! s; v  F7 p            d=d/a[v];
          }
5 h* V% X. j, P5 g) c" C6 U      }
2 W' U! B9 B! d8 m. ]/ t    return(2);
  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
  double x0,h,t,y[];
# [# t& ?2 r8 |  { int i,j,k,m;
  t  J% ~5 @. C, a    double z,s,xi,xj;
8 Q6 x: g& B# g( x$ N: C    float p,q;
    z=0.0;
+ H( X: w3 M' l    if (n&lt;1) return(z);
    if (n==1) { z=y[0]; return(z);}
# R. S$ D$ Z% F7 S  w" C    if (n==2)
. B& w: @! o$ `( b      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
$ _$ d6 B2 z) G  Z        return(z);
8 x+ r7 Q* p* i1 I      }
: k& _* |2 [, {2 z- t* P" Z& v5 U6 m    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
2 J) V0 q3 x) _" C9 E3 D2 s' o$ E3 a/ x        if (p&gt;q) i=i+1;
2 K3 p- I; ^. }" F* c0 P      }
    else i=0;
    k=i-4;
    if (k&lt;0) k=0;
    m=i+3;
" W* v$ K. ]" B) c! I    if (m&gt;n-1) m=n-1;
6 F' Q1 t3 E! n) a! X    for (i=k;i&lt;=m;i++)
8 f# R% I5 R8 i: E& @# P      { s=1.0; xi=x0+i*h;
5 K  o4 G. b  Z3 L/ J        for (j=k; j&lt;=m; j++)
          if (j!=i)
" K4 s3 V. F" Z9 q) n3 I9 ^/ K            { xj=x0+j*h;
              s=s*(t-xj)/(xi-xj);
# M  C5 q" {' b& ?! O: C' s            }
        z=z+s*y;
& {* a& y. K1 Q: t- `# g      }
    return(z);
; o* \) n6 @% w4 c& ]. i) ^; ^! {  }
2 f3 a$ n* a7 L向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
. A0 a( z) O9 [5 C2 p  int n,m;
  double x[],y[],a[],dt[];
; e! j4 q. N9 r# C, d; {$ z" x) ^2 X  { int i,j,k;
: P' `, |, j4 L5 E    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;
    if (m&gt;n) m=n;
$ q0 w+ t: a6 m4 t4 i$ s4 ^    if (m&gt;20) m=20;
    z=0.0;
# w  x3 b( @" t- S    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;}
4 C: _" }5 f# l* v6 O    c=c/d1; p=p/d1;
    a[0]=c*b[0];
( ~( N+ ^* L( e9 g0 S" ^* T    if (m&gt;1)
2 d4 S! M) A9 R- N' @% K, P0 l, R+ `      { t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
. ]8 @9 b; [* Q+ O- M- u. v          { q=x-z-p; d2=d2+q*q;
! A5 u8 K* w9 w1 \            c=c+y*q;
            g=g+(x-z)*q*q;
! d6 h0 h% J$ D5 l, M6 C! ^          }
6 \2 D% k8 ]  h1 v  K: b0 a5 p/ O5 f        c=c/d2; p=g/d2; q=d2/d1;
  B+ w) d0 |6 c6 ?. d        d1=d2;
        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
4 j% e6 E: R8 {/ k. l      }
    for (j=2; j&lt;=m-1; j++)
      { s[j]=t[j-1];
        s[j-1]=-p*t[j-1]+t[j-2];
        if (j&gt;=3)
; h' u3 q+ b0 o7 i+ q  `          for (k=j-2; k&gt;=1; k--)
2 [) P4 W6 g7 d0 {0 V            s[k]=-p*t[k]+t[k-1]-q*b[k];
        s[0]=-p*t[0]-q*b[0];
# M- }" l8 M# p; R        d2=0.0; c=0.0; g=0.0;
        for (i=0; i&lt;=n-1; i++)
6 i- @6 c. \) e* C6 C# r' k          { q=s[j];
            for (k=j-1; k&gt;=0; k--)
" N% v& o; Q3 l- U7 j4 @' w- l              q=q*(x-z)+s[k];
            d2=d2+q*q; c=c+y*q;
            g=g+(x-z)*q*q;
: t, h# `5 g$ M5 t. [) F/ v9 J          }
        c=c/d2; p=g/d2; q=d2/d1;
3 _# M4 U8 A! q: W  k        d1=d2;
        a[j]=c*s[j]; t[j]=s[j];
2 C0 r8 `. @- v# k4 U  T" Y        for (k=j-1; k&gt;=0; k--)
& ]% T$ ?6 N5 @5 c$ C9 g6 b+ I          { a[k]=c*s[k]+a[k];
" X/ V/ e  b" z3 i1 V. a            b[k]=t[k]; t[k]=s[k];
          }
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
, M& F6 f% N& H+ y! n) ~0 z- n2 [    for (i=0; i&lt;=n-1; i++)
+ N1 Z8 |. c# U9 [: w" y      { q=a[m-1];
' n1 X) N, K9 v( u7 a        for (k=m-2; k&gt;=0; k--)
          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;
4 p1 k/ A0 H7 O! I/ F' h8 h        dt[1]=dt[1]+fabs(p);
      }
1 q8 b9 t; s" l% c- `% g6 `    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
  double fromb(a,b,eps)
' a5 s" C3 K- }' W# ^3 x5 J% D  double a,b,eps;
9 F, R# I) t. c  { extern double frombf();
: G+ ^% ^( I+ w# S1 d2 I$ R    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
/ A' @8 R6 ^7 p! y9 A8 B( w9 o. |    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
    m=1; n=1; ep=eps+1.0;
$ R/ Y' l. [4 l7 C& f; C    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
/ X( c0 y* L$ I( f      { p=0.0;
        for (i=0;i&lt;=n-1;i++)
          { x=a+(i+0.5)*h;
            p=p+frombf(x);
5 D" J9 S4 |. w: N0 D4 u$ o3 s          }
# c- `0 N4 M, p! q/ M        p=(y[0]+h*p)/2.0;
        s=1.0;
/ u5 ~4 \0 E5 S# M& f' T        for (k=1;k&lt;=m;k++)
1 \& X' q# B5 ~3 T7 E$ p! Z7 [          { s=4.0*s;
4 V# u  x' Z! P/ ]- C            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
) e2 p" k! x. j2 Q& g1 L/ \5 t0 R0 B          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
/ {8 t' G- k% b3 C2 h      }
/ f' D( a! N( i5 j7 C5 |& N4 M& Z    return(q);
7 M  @7 [( `( r. N  }</P><>呵呵 希望对你有用！！</P>

---

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

---

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

---

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

Powered by Discuz! X2.5