数学建模社区-数学中国

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

---

作者: 西西    时间: 2004-5-10 08:46
标题: 求助！！ 数值分析程序（c语言）
<>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。</P>
<>1。牛顿迭代法</P>
5 _' e4 e# B; o<>1。全主消元法</P>
<>3。改进平方根法</P>
<>4。牛顿向前，向后插值法</P>
<>5。加权最小二乘法</P>
# X$ ]. Z* {  @9 |5 ~- a$ {0 t<>6。龙贝格积分法</P>
, R, W8 B* R" O) d' b% y$ M<>拜托了。。。。</P>

---

作者: loveyuan    时间: 2004-5-11 08:32
<>1。牛顿迭代法</P><>
9 m; z( Y" b5 Q  #include "stdio.h"
  #include "math.h"
5 v2 \( W6 F  f  int dnewt(x,eps,js)
( |/ G. U3 i# m6 T% G& z8 h  int js;
7 v/ ^. [3 |5 H+ q8 ]) K  double *x,eps;
3 w% h* v8 [4 B( Q& C2 i1 y  { extern void dnewtf();
    int k,l;
$ w* h2 s6 ^  g( t5 k+ \7 [    double y[2],d,p,x0,x1;
    l=js; k=1; x0=*x;
    dnewtf(x0,y);
    d=eps+1.0;
    while ((d&gt;=eps)&amp;&amp;(l!=0))
( W& L6 i: y7 F3 K6 e      { if (fabs(y[1])+1.0==1.0)
( T7 E' t; V0 X) U3 X          { printf("err\n"); return(-1);}
/ G. w7 k' A8 H6 ?# i        x1=x0-y[0]/y[1];
) f. }7 {# i4 c; P( Z        dnewtf(x1,y);
5 R8 X  O2 C. Z$ q& R        d=fabs(x1-x0); p=fabs(y[0]);
# @2 l/ q& Y2 I1 K. B- d2 Y- H        if (p&gt;d) d=p;
        x0=x1; l=l-1;
- F- m+ p9 _9 q# Z" `3 _8 i      }
+ d9 B% E: f6 I    *x=x1;
1 s4 K8 F8 \" U4 x2 x5 A$ t    k=js-l;
    return(k);
  }</P><>全主消元法</P><>#include "stdlib.h"
  #include "stdio.h"
  int acgas(ar,ai,n,br,bi)
5 m7 K# C. s: k) p  int n;
# r; K& z" M0 W9 y  double ar[],ai[],br[],bi[];
3 A; D2 v  Y$ o2 |  { int *js,l,k,i,j,is,u,v;
" {* `, \2 Q- ^$ T% }    double p,q,s,d;
; n% z/ P- _5 K    js=malloc(n*sizeof(int));
2 ?- f# f! i( {3 j6 e    for (k=0;k&lt;=n-2;k++)
$ e5 _. y* G2 z2 r) T5 Q) J      { d=0.0;
        for (i=k;i&lt;=n-1;i++)
# J! w% u( N' y; _        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;}
          }
7 }0 w* _; T4 v) V        if (d+1.0==1.0)
          { free(js); printf("err**fail\n");
# u) ]; R! t5 q) A6 v+ Q            return(0);
# T: w' u2 [: r) ^" A          }
        if (is!=k)
7 B/ M# o1 v. ]" s, I          { for (j=k;j&lt;=n-1;j++)
) h0 u3 P2 y) Z$ P              { u=k*n+j; v=is*n+j;
                p=ar; ar=ar[v]; ar[v]=p;
                p=ai; ai=ai[v]; ai[v]=p;
              }
            p=br[k]; br[k]=br[is]; br[is]=p;
- V' }0 t6 n1 v! A+ s1 S            p=bi[k]; bi[k]=bi[is]; bi[is]=p;
4 e, q; _/ R) M' p5 e          }
        if (js[k]!=k)
5 f+ i5 \7 b# ~3 T" G+ @/ }          for (i=0;i&lt;=n-1;i++)
            { u=i*n+k; v=i*n+js[k];
( M* R" t. J+ [. N              p=ar; ar=ar[v]; ar[v]=p;
% b+ I# T% d& n" O, ]) H              p=ai; ai=ai[v]; ai[v]=p;
            }
        v=k*n+k;
        for (j=k+1;j&lt;=n-1;j++)
; h& Q3 n5 v0 i! K4 A2 T! w" r          { u=k*n+j;
; N' R& T8 M" g& \            p=ar*ar[v]; q=-ai*ai[v];
, r" |* i9 h0 r            s=(ar[v]-ai[v])*(ar+ai);
; Z9 b; J; h5 h& A3 i" ?            ar=(p-q)/d; ai=(s-p-q)/d;
/ C, X* I9 D$ X# s* u" A! a          }
        p=br[k]*ar[v]; q=-bi[k]*ai[v];
7 |1 y- W9 ?2 r' v+ w5 R9 N* {3 x$ Z        s=(ar[v]-ai[v])*(br[k]+bi[k]);
' D' k4 l8 z9 q; {9 d6 \        br[k]=(p-q)/d; bi[k]=(s-p-q)/d;
        for (i=k+1;i&lt;=n-1;i++)
/ r$ ~5 J$ }; y          { u=i*n+k;
            for (j=k+1;j&lt;=n-1;j++)
              { v=k*n+j; l=i*n+j;
9 {" `; G: _) Y6 L/ g% s6 \, I" l0 `4 I                p=ar*ar[v]; q=ai*ai[v];
                s=(ar+ai)*(ar[v]+ai[v]);
  M3 G. A" X, H$ q' j7 U                ar[l]=ar[l]-p+q;
                ai[l]=ai[l]-s+p+q;
              }
7 U$ c6 x+ ~9 P1 X            p=ar*br[k]; q=ai*bi[k];
/ [$ Z7 \9 d. x% P$ a, b/ ]            s=(ar+ai)*(br[k]+bi[k]);
            br=br-p+q; bi=bi-s+p+q;
          }
/ f8 |, N: [* C  i1 @1 u      }
    u=(n-1)*n+n-1;
& R" p; W/ |% a4 m! x    d=ar*ar+ai*ai;
; a9 U0 P) c- A& P: J2 s/ w    if (d+1.0==1.0)
8 E+ [1 a1 {& K+ n/ p2 i* O* ~      { free(js); printf("err**fail\n");
( h  H8 `# Z3 x/ ]0 H" @        return(0);
      }
# b" x& [# [4 g    p=ar*br[n-1]; q=-ai*bi[n-1];
    s=(ar-ai)*(br[n-1]+bi[n-1]);
0 c1 A7 U1 u$ d2 h/ T: K    br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
    for (i=n-2;i&gt;=0;i--)
/ ~; A! O& U6 k) @" Y; J5 {    for (j=i+1;j&lt;=n-1;j++)
      { u=i*n+j;
# _, O  E5 X! ?; j# f1 b        p=ar*br[j]; q=ai*bi[j];
. y' @9 H& w* S( \+ |1 ]        s=(ar+ai)*(br[j]+bi[j]);
6 |: |' {7 X0 C" p' X& B) r8 u        br=br-p+q;
        bi=bi-s+p+q;
      }
1 d6 w* F2 b3 c  F5 ^    js[n-1]=n-1;
3 q: ?3 Q% u$ I    for (k=n-1;k&gt;=0;k--)
      if (js[k]!=k)
        { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
3 q, N  d- {# o7 ]5 k- U$ [( J          p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
        }
    free(js);
    return(1);
+ \/ I' ]& u( E$ ]2 H  }</P><>平方根法</P><>#include "math.h"
; f3 ?9 B3 @4 Q% g% I3 I  #include "stdio.h"
7 d( j# k$ b( o* `5 \  int achol(a,n,m,d)
  int n,m;
! b! X% s1 k) F3 @8 F: Z  p  double a[],d[];
  { int i,j,k,u,v;
) K/ h1 r+ S5 s( c; d    if ((a[0]+1.0==1.0)||(a[0]&lt;0.0))
      { printf("fail\n"); return(-2);}
    a[0]=sqrt(a[0]);
" g1 c9 r9 q; C8 K; j- [    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;
; D) W+ f- S: ?2 ~- K        for (j=1; j&lt;=i; j++)
          { v=(j-1)*n+i;
6 J$ e; O6 X) T3 f2 o9 f6 y. B            a=a-a[v]*a[v];
9 K7 z4 @. T2 O+ l8 b3 b; z5 b          }
5 z" E7 [3 H- z6 e" ~6 R" h% W        if ((a+1.0==1.0)||(a&lt;0.0))
          { printf("fail\n"); return(-2);}
        a=sqrt(a);
/ a5 G; h' _$ X5 ?        if (i!=(n-1))
& A/ M! f% s9 P# B2 _+ L+ V          { for (j=i+1; j&lt;=n-1; j++)
              { v=i*n+j;
8 ^1 H/ B2 E7 q9 n; x                for (k=1; k&lt;=i; k++)
  c- A/ r3 P# I6 H9 z9 O4 w" P7 ~                  a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
1 ^9 k0 ?! L1 q" h& w                a[v]=a[v]/a;
& k) L# K$ Q1 k% ^5 t              }
& t* O2 }5 |2 p+ m4 g2 {0 x          }
3 V, \/ c& _1 I3 q      }
    for (j=0; j&lt;=m-1; j++)
      { d[j]=d[j]/a[0];
, P* e% f* n$ c/ A/ j# F0 p0 F9 C5 h        for (i=1; i&lt;=n-1; i++)
          { u=i*n+i; v=i*m+j;
            for (k=1; k&lt;=i; k++)
0 `. d- X' w  F8 u              d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
1 E3 V. U9 Z: I" }3 v( G* G$ X            d[v]=d[v]/a;
          }
( z' h0 w& C' K8 D$ }      }
    for (j=0; j&lt;=m-1; j++)
      { 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++)
              { v=(k-1)*n+i;
                d=d-a[v]*d[i*m+j];
5 q! o  Z( v* R0 @4 J% g              }
* q8 J9 j2 ]% n; K8 w  c0 E/ ?            v=(k-1)*n+k-1;
            d=d/a[v];
          }
      }
$ R" P8 }& V8 B( v    return(2);
  }</P><>牛顿向前，向后插值法</P><>double eelgr(x0,h,n,y,t)
  int n;
6 n& F' w/ D+ C  double x0,h,t,y[];
+ s5 [1 n+ T, W, {6 n+ O2 a( P& b  { int i,j,k,m;
) Q9 J) v& i/ X" f    double z,s,xi,xj;
    float p,q;
    z=0.0;
    if (n&lt;1) return(z);
0 }  H. F$ Y* }4 m    if (n==1) { z=y[0]; return(z);}
    if (n==2)
; H: _3 t3 U5 U' N      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
) Q' u' E: S6 o        return(z);
      }
    if (t&gt;x0)
      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p&gt;q) i=i+1;
  r* E* K) K' g9 H" ~* w      }
( o2 t0 j0 \. I    else i=0;
! _4 t2 R" I6 X( `2 G. R( J$ {, o    k=i-4;
    if (k&lt;0) k=0;
# V2 N; X1 l. A- L# g. i    m=i+3;
) w" l; l# c8 i    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++)
          if (j!=i)
( G6 T( q: y; ~. E4 y            { xj=x0+j*h;
8 H. Q; E4 h6 Z) b6 ~              s=s*(t-xj)/(xi-xj);
% `% E, ~4 T' I4 O6 y7 j# u/ ^' @. g            }
' G5 h3 O' }( `5 ]        z=z+s*y;
9 x( Q5 `1 a. I. I( O      }
    return(z);
+ q+ L4 j$ {; C5 y: E7 l  }
向前，向后是一样的思想！！</P><>加权最小二乘法</P><>#include "math.h"
  void hpir1(x,y,n,a,m,dt)
2 ^* {. F/ W$ s( t$ y4 w; Z  int n,m;
% L0 m0 {/ G0 u: k# s8 d  double x[],y[],a[],dt[];
  { int i,j,k;
    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;
    if (m&gt;20) m=20;
    z=0.0;
2 B* Z2 Z: T* M0 \2 d% l. G    for (i=0; i&lt;=n-1; i++) z=z+x/(1.0*n);
3 s7 q4 v. t# `8 a  ^    b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
0 f  h6 F7 Z8 Y    for (i=0; i&lt;=n-1; i++)
# p# Q0 P6 n) G, V' z      { p=p+(x-z); c=c+y;}
    c=c/d1; p=p/d1;
    a[0]=c*b[0];
8 _) X! g* N6 O( D( h. [2 N    if (m&gt;1)
      { t[1]=1.0; t[0]=-p;
* [' d4 M; d3 |- ^5 F9 |        d2=0.0; c=0.0; g=0.0;
( h( ?# d2 {$ q; z        for (i=0; i&lt;=n-1; i++)
& H& f# x, X4 T# V* p) k          { q=x-z-p; d2=d2+q*q;
            c=c+y*q;
" x! o9 |' n8 S% c  D' B9 S# w1 ]            g=g+(x-z)*q*q;
' t4 p7 H; J! f5 [9 W1 n- P$ m          }
        c=c/d2; p=g/d2; q=d2/d1;
2 E5 B% ?! \" F/ A        d1=d2;
; H9 O8 x8 M; h! `$ k. }        a[1]=c*t[1]; a[0]=c*t[0]+a[0];
4 @( H; B. w: a9 F      }
    for (j=2; j&lt;=m-1; j++)
' v  q2 O3 V* h8 p$ {, j5 J      { s[j]=t[j-1];
        s[j-1]=-p*t[j-1]+t[j-2];
5 [- v1 |4 y! `3 O: f        if (j&gt;=3)
# U6 A. |* K# F- r          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];
8 @5 w# P0 k) [        d2=0.0; c=0.0; g=0.0;
* c7 M0 _8 b* P4 }* I        for (i=0; i&lt;=n-1; i++)
          { q=s[j];
            for (k=j-1; k&gt;=0; k--)
. h, V5 o/ e  X/ ]              q=q*(x-z)+s[k];
2 H( t' H; u2 S+ x0 ~& ]" N            d2=d2+q*q; c=c+y*q;
            g=g+(x-z)*q*q;
9 c5 v  Z+ U6 c/ N          }
        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];
          }
      }
    dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
" x0 t  C4 O) `$ ^    for (i=0; i&lt;=n-1; i++)
      { q=a[m-1];
        for (k=m-2; k&gt;=0; k--)
          q=a[k]+q*(x-z);
        p=q-y;
  U1 l$ ]6 v) O8 n2 L1 v- E        if (fabs(p)&gt;dt[2]) dt[2]=fabs(p);
        dt[0]=dt[0]+p*p;
        dt[1]=dt[1]+fabs(p);
      }
" H4 w& e7 y  c. P7 J  N6 A    return;
  }</P><>龙贝格积分法</P><>#include "math.h"
  double fromb(a,b,eps)
; v. M3 _) m' `8 ]$ m( Z8 M$ n  double a,b,eps;
  { extern double frombf();
2 T9 o% X. e( s# G+ W. K5 V4 M    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
    h=b-a;
7 ^: g9 e& `' E) e& j- |! V' K6 l    y[0]=h*(frombf(a)+frombf(b))/2.0;
$ J  n! X% |( {# d& w) a9 `6 u    m=1; n=1; ep=eps+1.0;
2 r# A7 R3 \* }+ |3 `  H    while ((ep&gt;=eps)&amp;&amp;(m&lt;=9))
; a7 `- ~* M4 S' ^) R/ O      { p=0.0;
  f" M0 q! V* L6 D5 [% {" X        for (i=0;i&lt;=n-1;i++)
          { x=a+(i+0.5)*h;
            p=p+frombf(x);
          }
        p=(y[0]+h*p)/2.0;
- T) M# O( v7 b4 s        s=1.0;
7 w' f6 b) H1 N4 T! q) k        for (k=1;k&lt;=m;k++)
. G" h0 Z. @/ N2 Q0 |          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
* {) A- r6 u" N0 h( l5 d% I6 Z            y[k-1]=p; p=q;
7 C) l" N# B& Q          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
3 G, q2 w! M" K% F, f6 t    return(q);
  }</P><>呵呵 希望对你有用！！</P>

---

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

---

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

---

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

Powered by Discuz! X2.5