查看: 6015|回复: 3

求助！！数值分析程序（c语言）

字体大小: 正常放大

西西

1 主题	0 听众	20 积分

升级 15.79%

该用户从未签到

电梯直达

1^#

发表于 2004-5-10 08:46 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

<

>本人紧急求助，希望各位大虾们伸出你们友情之手，帮我一下。
<

>1。牛顿迭代法
<

>1。全主消元法
<

>3。改进平方根法
<

>4。牛顿向前，向后插值法
<

>5。加权最小二乘法
<

>6。龙贝格积分法
<

>拜托了。。。。

zan

转播0 淘帖0 分享0 收藏0 支持0 反对0 微信

使用道具举报

loveyuan

21 主题	2 听众	104 积分

该用户从未签到

2^#

发表于 2004-5-11 08:32 |只看该作者 |招呼Ta 关注Ta

<

>1。牛顿迭代法<

>
 #include "stdio.h"
 #include "math.h"
 int dnewt(x,eps,js)
 int js;
 double *x,eps;
 { extern void dnewtf();
int k,l;
double y[2],d,p,x0,x1;
l=js; k=1; x0=*x;
dnewtf(x0,y);
d=eps+1.0;
while ((d>=eps)&&(l!=0))
 { if (fabs(y[1])+1.0==1.0)
 { printf("err\n"); return(-1);}
 x1=x0-y[0]/y[1];
 dnewtf(x1,y);
 d=fabs(x1-x0); p=fabs(y[0]);
 if (p>d) d=p;
 x0=x1; l=l-1;
 }
*x=x1;
k=js-l;
return(k);
 }<

>全主消元法<

>#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));
for (k=0;k<=n-2;k++)
 { d=0.0;
 for (i=k;i<=n-1;i++)
 for (j=k;j<=n-1;j++)
 { u=i*n+j;
 p=ar*ar+ai*ai;
 if (p>d) {d=p;js[k]=j;is=i;}
 }
 if (d+1.0==1.0)
 { free(js); printf("err**fail\n");
 return(0);
 }
 if (is!=k)
 { for (j=k;j<=n-1;j++)
 { 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;
 p=bi[k]; bi[k]=bi[is]; bi[is]=p;
 }
 if (js[k]!=k)
 for (i=0;i<=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;
 }
 v=k*n+k;
 for (j=k+1;j<=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;
 for (i=k+1;i<=n-1;i++)
 { u=i*n+k;
 for (j=k+1;j<=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;
 }
 p=ar*br[k]; q=ai*bi[k];
 s=(ar+ai)*(br[k]+bi[k]);
 br=br-p+q; bi=bi-s+p+q;0 F6 ~; v3 ^+ H
 }1 t( }* E) H3 {+ B# |
 }
( o& o' G; L$ L E% Q3 `: E u=(n-1)*n+n-1;/ e* P8 ?, c, |5 V$ u
d=ar*ar+ai*ai;
( a& I- H* r4 r8 i8 `1 Z# f( x3 l if (d+1.0==1.0)
. t7 v- s- D6 j { free(js); printf("err**fail\n");
$ P- \6 q6 i$ K+ m+ f5 d5 I# | return(0);
. U. M: Y2 T6 f ?4 V, s f }
4 e# o! {. j' b, z$ v p=ar*br[n-1]; q=-ai*bi[n-1];& a( K- K8 `8 x, e8 R! v8 _" H! B$ t
s=(ar-ai)*(br[n-1]+bi[n-1]);8 s/ d9 }/ ~* {* Y
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;: {" h+ p( J: t/ {2 q
for (i=n-2;i>=0;i--)
0 j, |( @+ |* t for (j=i+1;j<=n-1;j++). Q$ z. B/ q% s) M
 { u=i*n+j;9 \6 N8 u: j" I
 p=ar*br[j]; q=ai*bi[j];1 P; m3 Y/ V" ?* R3 j' `
 s=(ar+ai)*(br[j]+bi[j]);: |! s0 F/ f( o8 P) w
 br=br-p+q;
' h3 K3 q+ B8 ~. I8 W+ A7 O' y bi=bi-s+p+q;
 @( B3 D+ w" L/ U }
js[n-1]=n-1;
1 @+ e; `: ^! }5 c0 ~ for (k=n-1;k>=0;k--): Y, ~8 t2 m6 A. `. o1 [
 if (js[k]!=k)
8 H& d7 k8 k1 m6 k: O* I4 w, u1 U { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;9 I A% s5 ^# Z, I( Z& t
 p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
4 V2 {' _7 B! f5 J9 K! `! w6 G }
 ?% m: p0 K' m" D+ P1 D! B4 T) U free(js);
, m% h/ N U! P) y return(1);7 R2 u- ^9 G' u' ?
 }<>平方根法<>#include "math.h"
$ J) |+ D0 F W- R6 u #include "stdio.h"
; c/ E2 d% S; I" X6 H1 U8 q7 c1 q int achol(a,n,m,d)! t6 o' H- ]7 L) y
 int n,m;0 z" u' j2 f9 _# R |' b4 l/ r8 ^
 double a[],d[];
9 i$ e j! Y( ^6 U8 f' ?4 ^" z( L+ Z { int i,j,k,u,v;( ?- j% u. L6 t3 f. w
if ((a[0]+1.0==1.0)||(a[0]<0.0))
9 ?- _* c1 `- g) {+ h { printf("fail\n"); return(-2);}/ z$ E Y# B4 p/ W/ {. D# d' Q
a[0]=sqrt(a[0]);
( b4 V: {7 M0 Z5 I for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];( L: N# h% K1 o
for (i=1; i<=n-1; i++)1 o+ f. T. k1 @$ e1 Q
 { u=i*n+i;
, F3 y$ q- [7 X2 r* [% G for (j=1; j<=i; j++)/ K# [* v: d/ o; G4 e) z
 { v=(j-1)*n+i; W; m9 {# p7 X6 `, n, J) e( s
 a=a-a[v]*a[v];
1 n& S) d5 |6 O; T }
9 Q7 ^$ l4 X9 @& d0 Q9 D2 x$ h1 ?6 U if ((a+1.0==1.0)||(a<0.0))6 s) ^' I6 r$ X8 U
 { printf("fail\n"); return(-2);}
1 V I( @+ g ~5 k" P1 z& @ a=sqrt(a);
0 H% F' o) n8 `5 N' Q- P4 \! q9 ~ if (i!=(n-1))& o" K, H* o' r1 M h
 { for (j=i+1; j<=n-1; j++)
, D; H5 Y1 T4 X8 a3 L$ G9 H* w { v=i*n+j;
- E7 c" m) E) `0 q7 X for (k=1; k<=i; k++)
1 ?. S/ {4 p7 ~) D a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];: t2 S/ x; b3 D) H3 y) V
 a[v]=a[v]/a;
8 D! G. k) e& y- n: ?" L }, _: V& n: i0 t9 A J" @
 }& I2 u8 ~: C8 k0 a" H
 }' i- k7 U, U$ }3 o; `# n% ~. f
for (j=0; j<=m-1; j++)" b7 Z& N2 f j
 { d[j]=d[j]/a[0];
 _6 H. g3 k. M9 Z" R0 U" O for (i=1; i<=n-1; i++)3 E2 u+ z2 I; I, s
 { u=i*n+i; v=i*m+j;/ L1 T, B" _3 Y) S. g V( @
 for (k=1; k<=i; k++). C, D2 E& L0 Q* o. c
 d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];; P+ T0 ~, _5 k' ^/ l( p i
 d[v]=d[v]/a;
$ {& d+ q+ J% p U' J }
2 ]( Q9 T3 Q+ \ }
: I- X4 `# B- }9 g: s for (j=0; j<=m-1; j++)
- y/ a/ }+ i# L! c$ g" c8 _: E { u=(n-1)*m+j;
. l$ g6 t% \' F9 a8 a d=d/a[n*n-1];! u2 K/ l& E9 r, T% h" e
 for (k=n-1; k>=1; k--)
" v) r9 W4 A+ U+ [ { u=(k-1)*m+j;
 p, C$ P( `% |4 e+ R" L for (i=k; i<=n-1; i++)! a4 d7 p. X& X" y8 ?8 m4 M$ h
 { v=(k-1)*n+i;
6 W% b' T r- P* p d=d-a[v]*d[i*m+j];8 k1 [2 S% |7 H) }2 X4 k
 }6 y N( D5 r* F
 v=(k-1)*n+k-1;3 {2 k5 g; `0 n4 H
 d=d/a[v];
4 Y, b& m9 `3 c/ Q% Q }
4 F8 Q3 a7 |9 I) P9 |# A1 z4 ` }/ B" }$ M) k: ]$ D( n' ?! T z
return(2);' n+ W! T) S8 {; B
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)2 g8 F: n. ?3 y( s) Z i3 b
 int n;
$ }% Y6 P2 n/ ]5 v# i double x0,h,t,y[];
# m8 }# W; Y/ ^ { int i,j,k,m;
( c. e' E/ [$ X, u+ _ double z,s,xi,xj;) L/ G2 }* S" v$ g' C+ s! ~
float p,q;
- ]6 z$ l7 N! X+ G1 d' _ z=0.0;) a% U p- P3 k0 Q. t6 q$ u; j
if (n<1) return(z);( p. q/ H/ ]" a! W) l! `/ d
if (n==1) { z=y[0]; return(z);}! [, v; N: d+ ~9 \
if (n==2)& c w: a6 h0 I/ N- r
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
. [1 V4 g1 h' t5 u5 k return(z);
: d4 [8 b( k3 ~% s" J) @( q }1 e* B+ ?; q! W3 _
if (t>x0)& z6 A0 k# r4 U, |
 { p=(t-x0)/h; i=(int)p; q=(float)i;
0 T# o" G7 \0 y7 z% u& e if (p>q) i=i+1;
& K7 {$ C" _: m* g }( z- I+ y6 m- m* T8 T5 Y Y e: A
else i=0;4 u. K0 B. m6 r! z: N* z
k=i-4;
+ N U' ?0 M# P: m6 X& Y& J; F if (k<0) k=0;$ v3 |' w) b* n H3 u, B7 r6 N4 C
m=i+3;
4 L. s, u4 x: e) M if (m>n-1) m=n-1;+ S0 E' g2 l- q$ i/ t x( _
for (i=k;i<=m;i++)9 o; u; l: m$ w* T2 T
 { s=1.0; xi=x0+i*h;
2 r. H% l6 O# i* w- X* U3 J+ o8 R; { for (j=k; j<=m; j++)
+ Z+ e5 x/ \& G! y/ W9 C if (j!=i)4 ^0 x4 q3 M1 o! j% t
 { xj=x0+j*h;
, n, M3 B9 C: x$ Y s=s*(t-xj)/(xi-xj);- v( [* x# a0 H( r- f* y
 }. P6 m, T, `9 k* c# I$ q3 Z8 e* D
 z=z+s*y;" w) d( L3 }- F
 }
. v- R1 W. i. n0 K return(z);7 B3 Y) I8 u: p" ]
 }
& P, \! C) b' V) A+ G& N向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"0 r4 W7 f" H+ ~6 |
 void hpir1(x,y,n,a,m,dt)% a1 x: M- G& z6 s) A2 B4 L1 q, K
 int n,m;
: ^) J0 l- `* x9 [! Y6 G4 |7 K double x[],y[],a[],dt[];
( R" R- e: F3 Q. b3 f3 y { int i,j,k;
# N& O, T1 W0 r) P" w3 h1 y double z,p,c,g,q,d1,d2,s[20],t[20],b[20];+ b2 g1 a" o0 y8 c4 f2 B, j
for (i=0; i<=m-1; i++) a=0.0;* |! Q9 O( m1 M# G% f. k
if (m>n) m=n;
6 I/ H' |/ _* O7 n3 L8 p5 J/ N if (m>20) m=20;) L n$ O2 D' A O
z=0.0;
" [8 b4 u- ?3 J2 ?: z/ o( Z for (i=0; i<=n-1; i++) z=z+x/(1.0*n);, t0 T4 x2 G, O$ `9 [
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
0 {0 ]* [8 ]5 T7 N" Y9 q8 k for (i=0; i<=n-1; i++)' e( B& b: u' B0 i
 { p=p+(x-z); c=c+y;}1 o1 l2 w! s" K" Q8 d8 z% _
c=c/d1; p=p/d1;
* S/ O3 f% _1 \ a[0]=c*b[0];( S% ~3 E! [8 Z6 Q( g$ @. \& N
if (m>1); t4 ^; }: b8 L
 { t[1]=1.0; t[0]=-p;
# l p1 [4 T4 t% k d2=0.0; c=0.0; g=0.0;
, C0 b4 q4 a# y3 z7 k for (i=0; i<=n-1; i++)
" b2 t4 i6 m/ a { q=x-z-p; d2=d2+q*q;+ c0 j+ _8 D$ }8 {
 c=c+y*q;$ l/ |( K6 s5 E- X+ }/ J
 g=g+(x-z)*q*q;! x" r& ?8 p5 d" j. }& j: ^9 e5 J
 }
4 \* k- X, |+ i6 G# {4 w! L( k c=c/d2; p=g/d2; q=d2/d1;7 f: H% i4 c' h! w0 C; e
 d1=d2;
4 a% I5 h+ G( Q a[1]=c*t[1]; a[0]=c*t[0]+a[0];
2 a' K: G2 p7 k }5 Y) R. S) s0 O( G: g9 O4 P
for (j=2; j<=m-1; j++); `) J3 {: N9 i0 e/ }7 _
 { s[j]=t[j-1];9 Q8 K4 R3 p' z
 s[j-1]=-p*t[j-1]+t[j-2];" y1 N* `- u" n+ b4 _# j; z' h
 if (j>=3)5 c a& Q' d" V. K
 for (k=j-2; k>=1; k--)& _5 `% x" T! e+ C# `
 s[k]=-p*t[k]+t[k-1]-q*b[k];
* K4 a* `' @" ]+ Z4 v s[0]=-p*t[0]-q*b[0]; l$ s& S! V4 ^* ~& ~
 d2=0.0; c=0.0; g=0.0;4 n6 M5 P) s {+ m3 b
 for (i=0; i<=n-1; i++)
& | ?0 X/ i# M8 D! { { q=s[j];5 ?3 B4 U! U0 d& |. B0 k, Z
 for (k=j-1; k>=0; k--)( d( P# w ~# ?% N+ C
 q=q*(x-z)+s[k];! s5 A4 b* d+ R: a
 d2=d2+q*q; c=c+y*q;8 o( F% Y3 T5 E2 g# \0 E* B
 g=g+(x-z)*q*q;
 g: i& v: J; Z+ s }
, k7 f3 C3 V- l# d! [ c=c/d2; p=g/d2; q=d2/d1;
/ g* E8 y0 d0 n& ~1 s# K d1=d2;* m; Z- X1 d( g, \ F D2 Z- i
 a[j]=c*s[j]; t[j]=s[j];
# @3 {8 p- v) ~* J+ o+ S* ]# a+ J for (k=j-1; k>=0; k--)
% I ^" ~1 x6 s, N3 L7 q { a[k]=c*s[k]+a[k];# | j0 y& j7 k1 K" p* d% j
 b[k]=t[k]; t[k]=s[k];
3 w: [" X; c A6 H6 i$ U$ J }
" B% C; H- T$ i5 `5 @3 T0 u' i% Y }
/ f/ l) E4 Z0 Y8 |& N dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
" V' X3 a# X1 w/ w, o( W' C for (i=0; i<=n-1; i++)2 \' y% q( r; _. J: z- n6 K
 { q=a[m-1];
- e/ ^ {& p9 {( b for (k=m-2; k>=0; k--)2 K* Q) {1 C+ j: C
 q=a[k]+q*(x-z);; Y( C1 n- y# z4 [, I& p* \) t m
 p=q-y;
' {6 ?, s6 p/ m if (fabs(p)>dt[2]) dt[2]=fabs(p);6 e' l* {% Q% u1 b3 D6 U! I8 C
 dt[0]=dt[0]+p*p;* Z. U0 j9 _% E+ g# E
 dt[1]=dt[1]+fabs(p);
2 |: Q7 A1 k/ [$ D }( T" T" R" U4 u) `0 v6 U
return;, }' x0 C* P, b. s& t& ?' v- H9 ?
 }<>龙贝格积分法<>#include "math.h"+ L/ ^7 }$ \* V( u* k5 M
 double fromb(a,b,eps), Q0 U4 C' W& k& ^0 x
 double a,b,eps;
4 y# J" F! w Q3 V1 [+ |2 Z8 \4 D { extern double frombf();7 X! J' `. T; s- c/ q# |/ ?& k
int m,n,i,k;3 ^% M- W. H6 [% n# N2 y3 |) ]
double y[10],h,ep,p,x,s,q;* L+ I! o, ?+ Z( A5 r% e n
h=b-a;
# A0 L& U2 K) B y[0]=h*(frombf(a)+frombf(b))/2.0;' a3 W! n0 x' c3 z9 H. }; q9 L
m=1; n=1; ep=eps+1.0;0 I! f, _# P2 k! ~
while ((ep>=eps)&&(m<=9))# `) T9 q& u0 E1 m
 { p=0.0;
* w# j5 [* @7 j/ n: F. E4 Z! Y for (i=0;i<=n-1;i++)
6 G% W% M' F% V/ ]! |+ \( D { x=a+(i+0.5)*h;( x, P2 e+ s) o- R; `; | T$ C
 p=p+frombf(x);
! |) v$ P* f- [2 G. r: N }1 _0 Z4 a7 L/ r% T9 D
 p=(y[0]+h*p)/2.0;
: ~# ~5 Y$ [2 K7 {9 _ s=1.0;
: y/ W. M* {% ]0 }4 ~4 h# \ for (k=1;k<=m;k++)
9 t! G0 L+ V; w ], E/ z4 l. L { s=4.0*s;
0 S# g% G& j/ c q=(s*p-y[k-1])/(s-1.0);; L+ [; T, ]! O3 Y5 r
 y[k-1]=p; p=q;
. n" }/ z* V( ? }
1 H/ H( D, w& F ep=fabs(q-y[m-1]);
" e" C! g* n( P' U: R m=m+1; y[m-1]=q; n=n+n; h=h/2.0;. m$ f; R/ O' Q! \& r
 }1 {( M$ r+ e @3 _7 _$ x+ H+ r
return(q);( Z% c; d: C5 E: @* L) c
 }<>呵呵希望对你有用！！