数学建模社区-数学中国

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

---

作者: 西西    时间: 2004-5-10 08:46
标题: 求助！！ 数值分析程序（c语言）
<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
<>1。牛顿迭代法</P>
: ?! L# ^6 Q; V9 P( d<>1。全主消元法</P>
2 u: |2 v: r9 m- A( s<>3。改进平方根法</P>
5 Z5 E6 U9 C$ B; C' R0 a" h<>4。牛顿向前，向后插值法</P>
1 y! i- N4 ]/ d, s* \<>5。加权最小二乘法</P>
<>6。龙贝格积分法</P>
<>拜托了。。。。</P>

---

作者: loveyuan    时间: 2004-5-11 08:32
<>1。牛顿迭代法</P><>
  #include "stdio.h"
  #include "math.h"
1 E0 T9 @6 n7 I$ n: a% z# T  int dnewt(x,eps,js)
$ [$ a0 [8 N/ Z  Z  int js;
% Y' A5 b6 Z: |: z8 [) s  double *x,eps;
  { extern void dnewtf();
% F; {# s- E; J0 F    int k,l;
    double y[2],d,p,x0,x1;
) X9 ^" _6 p. o    l=js; k=1; x0=*x;
3 g6 B) R& @% B0 w    dnewtf(x0,y);
. ?$ y  V5 o8 u4 G# b    d=eps+1.0;
+ ~4 W& t# x1 O8 x# {6 n7 ?2 q& P% y: F0 ]    while ((d&gt;=eps)&amp;&amp;(l!=0))
      { if (fabs(y[1])+1.0==1.0)
          { printf("err\n"); return(-1);}
1 n) D: i4 C: b- \        x1=x0-y[0]/y[1];
1 U6 y, G+ M  J9 f8 G5 w4 D9 }5 X5 L        dnewtf(x1,y);
        d=fabs(x1-x0); p=fabs(y[0]);
        if (p&gt;d) d=p;
        x0=x1; l=l-1;
/ d& x7 r" v/ e4 k/ U2 [      }
5 O5 m: ]' P( @+ R    *x=x1;
5 D6 ?) J5 \: \1 B( y; _% ]    k=js-l;
2 @- x0 {/ w6 U6 M# A) D. j    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
& ~; z, S+ T5 Q2 S& X  int n;
: T2 Q# r# D- L$ I* t# K) r  double ar[],ai[],br[],bi[];
  { int *js,l,k,i,j,is,u,v;
; V! I. {8 @# F, J- k    double p,q,s,d;
& f# \9 h* p1 [    js=malloc(n*sizeof(int));
2 r) f* g* h7 L/ R8 d- a6 G    for (k=0;k&lt;=n-2;k++)
      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
+ m2 U) w. k0 Q. ?. c        for (j=k;j&lt;=n-1;j++)
9 @. i2 j+ }5 U6 j9 Q6 `0 Y' n          { u=i*n+j;
* J, ?; W' ?1 R1 g/ q            p=ar*ar+ai*ai;
            if (p&gt;d) {d=p;js[k]=j;is=i;}
6 a/ w, p8 Z9 r) T5 h          }
        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
* `( a: t3 m* C            return(0);
          }
        if (is!=k)
- x0 c- O* \2 s* Y' _& Z, K0 E2 k          { for (j=k;j&lt;=n-1;j++)
              { u=k*n+j; v=is*n+j;
" o% \1 e1 A+ A  Z                p=ar; ar=ar[v]; ar[v]=p;
- G1 y. e& I8 F# \$ u( g( S- y                p=ai; ai=ai[v]; ai[v]=p;
              }
            p=br[k]; br[k]=br[is]; br[is]=p;
            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
          }
3 `( x6 }" Z8 `8 `, ~        if (js[k]!=k)
          for (i=0;i&lt;=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
% i5 t& c, d7 e4 m8 j, \              p=ar; ar=ar[v]; ar[v]=p;
/ a! W; [4 S0 f! L! B; S0 p              p=ai; ai=ai[v]; ai[v]=p;
            }
        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
          { u=k*n+j;
            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]);
        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
* C) ?1 [# Q8 |* L& ^* H        for (i=k+1;i&lt;=n-1;i++)
: ~7 y- j: [3 I$ R/ L          { u=i*n+k;
) k: ]; @: \- B- d, F8 \            for (j=k+1;j&lt;=n-1;j++)
% B/ _0 Q! E9 Z4 K7 U              { v=k*n+j; l=i*n+j;
/ \" d1 f% l' Y7 L1 _( g; W                p=ar*ar[v]; q=ai*ai[v];
  F) f" {, q6 Z. a  r" g                s=(ar+ai)*(ar[v]+ai[v]);
                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
  c% w: H. n. L( W9 Q8 e            p=ar*br[k]; q=ai*bi[k];
            s=(ar+ai)*(br[k]+bi[k]);
5 t% P  i% N* b; l. w6 h4 j            br=br-p+q; bi=bi-s+p+q;
          }
      }
7 `( q9 e# Y, o* D# m4 j7 [    u=(n-1)*n+n-1;
4 s0 Y( ^$ k4 M$ ^    d=ar*ar+ai*ai;
5 n- W; P. v9 `2 A) |: e; e    if (d+1.0==1.0)
      { free(js); printf("err**fail\n");
        return(0);
      }
    p=ar*br[n-1]; q=-ai*bi[n-1];
9 r2 c; y5 Y! [: l' L7 e' z    s=(ar-ai)*(br[n-1]+bi[n-1]);
- {' J5 B" j3 X" D1 t" N    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
    for (i=n-2;i&gt;=0;i--)
- N' z- Z; F! Z, L# H0 m% g    for (j=i+1;j&lt;=n-1;j++)
% x5 n6 F! P* y2 P' O      { u=i*n+j;
* g2 {3 _& k& X4 P- y        p=ar*br[j]; q=ai*bi[j];
% g& C! `' J4 ~        s=(ar+ai)*(br[j]+bi[j]);
' ?1 d$ C9 Q. R" d7 Y        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;
0 I2 N; |  n* U! N$ h' j& h          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
        }
    free(js);
    return(1);
  }</P><>平方根法</P><>#include "math.h"
  #include "stdio.h"
8 u, r. o6 ]- D' V; _  int achol(a,n,m,d)
  int n,m;
; R0 Q# o2 S. r2 j) y  double a[],d[];
% e6 W( t4 A5 ]& @6 q  o6 Z$ k  { int i,j,k,u,v;
! W# o! M& g7 r7 a' ?  y/ e+ b    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
. m# s9 r, _8 N7 [4 j      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
& Y" `7 ]7 G' [! A. n    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;
% \* u5 t% l% ^: B. P9 ?# G' E        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
            a=a-a[v]*a[v];
          }
" _) U7 l0 o7 Y: f; A$ ^. i, w        if ((a+1.0==1.0)||(a&lt;0.0))
! M. }& H5 {2 X% G+ j          { 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;
; ^# D9 d6 k5 S& M5 @                for (k=1; k&lt;=i; k++)
4 r- Y  B+ T$ N  ?2 }& i                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
' w# g2 R! }4 H                a[v]=a[v]/a;
              }
          }
      }
    for (j=0; j&lt;=m-1; j++)
( ?: G/ k, N1 q. f      { d[j]=d[j]/a[0];
        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
            for (k=1; k&lt;=i; k++)
              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
4 I' |! ]' [3 X1 z6 R8 ^1 ?$ r            d[v]=d[v]/a;
          }
( B' T' L, u7 \0 F! ^      }
    for (j=0; j&lt;=m-1; j++)
! g7 U& ^" x& E$ O( m( b      { u=(n-1)*m+j;
; [; i+ M$ N2 S; t        d=d/a[n*n-1];
' w3 p1 z; E% _8 X        for (k=n-1; k&gt;=1; k--)
          { u=(k-1)*m+j;
7 C! ?. i; n. T  }; _8 t            for (i=k; i&lt;=n-1; i++)
8 I8 Z' V9 g1 p, z, q              { v=(k-1)*n+i;
                d=d-a[v]*d[i*m+j];
/ {# ^* d, I* B# ?. Q9 X              }
9 k- G4 J2 N1 c! J* r/ E1 i            v=(k-1)*n+k-1;
            d=d/a[v];
1 V: @3 }! V9 T* y          }
      }
    return(2);
% E- s4 i0 v- v1 [8 ?  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
# v4 |. ?4 t8 i  int n;
5 I7 T' g( ^5 A: m1 q  double x0,h,t,y[];
5 m2 S  a  H: f. N/ ~! Y" L  { int i,j,k,m;
    double z,s,xi,xj;
    float p,q;
4 v7 e0 w$ X3 i2 _, n/ Y! d  _    z=0.0;
    if (n&lt;1) return(z);
1 D6 \' k- ?( m' u0 r0 h    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return(z);
      }
: C0 H: X% A4 ]# a1 ?% [    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
  y: c5 Q& @: m- [  g7 Q+ S      }
    else i=0;
    k=i-4;
" S  h2 `' F* K2 A    if (k&lt;0) k=0;
: i# T7 {! A% D    m=i+3;
! x$ N- ?7 a6 `6 |# f- `    if (m&gt;n-1) m=n-1;
    for (i=k;i&lt;=m;i++)
      { s=1.0; xi=x0+i*h;
3 z. P8 P3 I/ o3 p+ m% g        for (j=k; j&lt;=m; j++)
0 K: E/ w# x6 h8 m; N          if (j!=i)
. m) l$ x5 Y7 J: L' M, |# |. M            { xj=x0+j*h;
8 B+ V$ K% |$ B  a3 A# T; A              s=s*(t-xj)/(xi-xj);
5 `) K2 T; b! P0 F            }
        z=z+s*y;
$ P2 ^, K! c1 d  G+ X8 p8 y5 d1 @      }
    return(z);
+ l+ K/ S! U0 e0 j" g2 x  }
& e( l6 o) R8 u向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
& H+ U) e9 @& _6 y2 g( ~5 }  int n,m;
$ u8 M' {$ p* ?- i4 Q7 B( ~4 |  double x[],y[],a[],dt[];
! V6 b8 R2 i$ l  { int i,j,k;
    double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
% F% r6 i( S* D$ ^    for (i=0; i&lt;=m-1; i++) a=0.0;
3 d) G6 u) ^" b$ _+ b" H+ R' Z    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);
; U3 X4 ?! k2 x8 L    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;
    a[0]=c*b[0];
) g6 V% S8 z# ?. R7 K$ g' Y    if (m&gt;1)
& t: f+ c  ?* T+ a      { t[1]=1.0; t[0]=-p;
        d2=0.0; c=0.0; g=0.0;
9 M3 e/ v! W+ `; @        for (i=0; i&lt;=n-1; i++)
. H' j7 T, B) C3 O          { q=x-z-p; d2=d2+q*q;
3 m0 j$ _! |+ N% e- t            c=c+y*q;
  C) [. b5 E1 ~8 `            g=g+(x-z)*q*q;
! l% s* P! [& K( g$ ]          }
        c=c/d2; p=g/d2; q=d2/d1;
% T! H0 \; D7 n- a" t        d1=d2;
        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
      }
& k7 b8 X8 n! W$ X+ `2 ^$ u% a  y    for (j=2; j&lt;=m-1; j++)
      { s[j]=t[j-1];
        s[j-1]=-p*t[j-1]+t[j-2];
" k5 W9 Z1 c/ E- M        if (j&gt;=3)
3 W! F2 B  T. Y. ~8 i# D          for (k=j-2; k&gt;=1; k--)
            s[k]=-p*t[k]+t[k-1]-q*b[k];
' g$ o. f" {: g* U+ o' }        s[0]=-p*t[0]-q*b[0];
: E# ^. T5 k, Q4 J3 v0 f        d2=0.0; c=0.0; g=0.0;
; V( ?3 }7 \) ]8 i. x- F        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
5 y& Z! |: _: l            for (k=j-1; k&gt;=0; k--)
. i- G9 t+ C1 a              q=q*(x-z)+s[k];
            d2=d2+q*q; c=c+y*q;
0 d7 [- Y. p) M            g=g+(x-z)*q*q;
; C4 d( |5 d1 j! G          }
2 c# L: g5 A- o        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--)
8 Q. S8 S# U+ c          { a[k]=c*s[k]+a[k];
            b[k]=t[k]; t[k]=s[k];
* E8 ^5 b7 w; N* M5 P2 L: f# u: \          }
. Z  H3 @# a7 c1 Z. {. m1 \      }
; `# |1 v% }: \    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
' Y) f2 S1 s1 |7 V- n    for (i=0; i&lt;=n-1; i++)
      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
0 E+ c9 k% P- ~( i  |& |          q=a[k]+q*(x-z);
5 R$ f. Q4 b  i( L        p=q-y;
        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
6 F. |; ~2 M5 P+ u& |        dt[1]=dt[1]+fabs(p);
* @$ h) H* P3 P8 U. j9 ]      }
    return;
- @3 Q! [1 `; w  B# {  }</P><>龙贝格积分法</P><>#include "math.h"
  double fromb(a,b,eps)
7 W; @) s$ u% o# [) G6 D! p; [  double a,b,eps;
+ g+ |' k1 ]) _; n7 t  { extern double frombf();
( r* l7 \, I$ K# g    int m,n,i,k;
: Z' z0 N# m% [, I2 o    double y[10],h,ep,p,x,s,q;
: [( t, o0 h( p: k    h=b-a;
    y[0]=h*(frombf(a)+frombf(b))/2.0;
) N8 ]! j7 M2 A: s% U0 @    m=1; n=1; ep=eps+1.0;
    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
      { p=0.0;
; A3 e6 o) ?" e/ K5 b7 x        for (i=0;i&lt;=n-1;i++)
          { x=a+(i+0.5)*h;
- `5 ?2 h5 I) d- o9 \            p=p+frombf(x);
          }
. M8 \+ G% x; L1 u        p=(y[0]+h*p)/2.0;
9 C& e' Q- _5 V. x- b0 Y        s=1.0;
7 o1 i; y9 a5 _, M        for (k=1;k&lt;=m;k++)
          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
( d2 d! q. h! q            y[k-1]=p; p=q;
          }
        ep=fabs(q-y[m-1]);
4 h$ D" O. u2 A! _9 _7 G. }" J        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
    return(q);
  }</P><>呵呵 希望对你有用！！</P>

---

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

---

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

---

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

Powered by Discuz! X2.5