查看: 6020|回复: 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;* j+ K# L' Y+ g& d0 M
 }- F# i/ n; d& }. z9 A) T
 }& Q8 `% N+ d5 e. |, {2 d
u=(n-1)*n+n-1;
; X$ t9 @+ J3 f3 p d=ar*ar+ai*ai;6 y) F+ t* ?* X& Y, C2 h( j
if (d+1.0==1.0)+ b2 }: ]! ~3 H7 T1 K: [4 X
 { free(js); printf("err**fail\n");6 S4 q, L$ c* L7 t6 |+ V( a$ Y
 return(0);
, m- n, A* W7 P. z }1 g! |, E e2 s" D/ A% ?, ]& c
p=ar*br[n-1]; q=-ai*bi[n-1];
s=(ar-ai)*(br[n-1]+bi[n-1]);
. k5 @# Y, F/ F3 ~% Q6 @2 c& }! | br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
" B# J7 _/ n/ k: C) @7 `4 D, B for (i=n-2;i>=0;i--)3 G7 y8 ^) u6 ^: |
for (j=i+1;j<=n-1;j++)& R% L) K3 l: D8 r! p
 { u=i*n+j;6 p0 `' G8 g6 e4 f
 p=ar*br[j]; q=ai*bi[j];* P" Q$ p% W: z4 [2 y
 s=(ar+ai)*(br[j]+bi[j]);5 l. o w, p+ \8 O
 br=br-p+q;5 @. Y. m) f1 b$ q+ I
 bi=bi-s+p+q;
8 m+ m2 H% t! ~& O }
$ y/ Z, R* M! {9 V& o9 T1 ]1 Q js[n-1]=n-1;5 e, | e. ^* y1 d/ ?0 N8 q0 G9 _
for (k=n-1;k>=0;k--)% g2 _4 j5 m7 d& I% k! s
 if (js[k]!=k)0 E2 ^, x3 ~- a5 }: r
 { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
( z$ U4 @! j& |2 J6 q p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
& Q4 o( G& l; k* s }
* W$ }8 D5 O7 T9 h' { free(js);
7 e( n! L9 z. M8 z return(1);2 Y% M$ I% d& O+ u5 h8 s+ P
 }<>平方根法<>#include "math.h"
7 D0 @* t) g6 L5 ?5 F) l; t #include "stdio.h"8 k( [2 u4 {1 m
 int achol(a,n,m,d)
 G% p \- X- T7 O! K6 ^8 k int n,m;
# ?0 e' W y3 {7 L* ^3 q: H2 E double a[],d[];8 Z1 q; `% G1 d& q
 { int i,j,k,u,v;
4 L: `3 Z+ @2 ~4 Y; j( C if ((a[0]+1.0==1.0)||(a[0]<0.0))% {; E0 G# g$ N5 G1 F
 { printf("fail\n"); return(-2);}6 U i( z6 |; B: n, U# C
a[0]=sqrt(a[0]);3 ^2 u5 B' l7 U Z, ]
for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
. m9 J% m+ `" ]# @+ j for (i=1; i<=n-1; i++)8 W$ x! j6 U# ~0 {3 E# [/ a
 { u=i*n+i;
0 H+ |+ @0 ]1 J# R- ~5 x/ o& h! @ for (j=1; j<=i; j++)
- A, z. ?* ~# |) d4 g { v=(j-1)*n+i;8 M3 E8 [3 t5 E% V& f/ h ?6 I$ h3 i
 a=a-a[v]*a[v];
% R5 M* ]: G* d }
6 E) x. f6 n& S if ((a+1.0==1.0)||(a<0.0))
% ~' R9 e- ] R" u { printf("fail\n"); return(-2);}
7 i" W( Y# v/ \: ~ a=sqrt(a);
9 l9 s9 @, \6 _+ m if (i!=(n-1))* r$ z. o& m- l4 f" ^
 { for (j=i+1; j<=n-1; j++)0 P9 |) ]0 v2 f
 { v=i*n+j;9 x1 G1 O+ d! @' b! m2 M3 B# g
 for (k=1; k<=i; k++)2 l5 `# D3 U& ` s/ {# h
 a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
3 ]. }! j/ H: f( @2 ]0 `; \ a[v]=a[v]/a;& G/ Y- R6 u& h- m4 T
 }
5 [, ^) \5 h) ~* U' g5 f }
3 [4 |( f U2 |9 T }
6 \4 M; D7 w& ^8 f6 ~9 { for (j=0; j<=m-1; j++)
" l, }' q3 h3 K6 l. v# i- } { d[j]=d[j]/a[0];7 s! F% ^, i3 C8 w Q6 g# Q/ o U
 for (i=1; i<=n-1; i++)7 P2 P6 X' ? \- l& W
 { u=i*n+i; v=i*m+j;
2 Z" v/ R( n, W for (k=1; k<=i; k++)0 x5 S6 D1 J1 W& I
 d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];; w$ l; }, r3 I% [* E+ _
 d[v]=d[v]/a;0 g+ L/ g4 N. Q" V$ X: E
 }3 f8 I$ o4 N. I: |9 E5 o; [- G
 }
% r. t7 ]" u9 L/ v for (j=0; j<=m-1; j++)
# p2 S' w# ?3 a* @' R# ~0 W { u=(n-1)*m+j;$ Q& ~9 f7 o! r0 `9 k
 d=d/a[n*n-1];
$ [: k4 o- q) C, n5 i; h5 r for (k=n-1; k>=1; k--)
) X ~; y( f$ x { u=(k-1)*m+j;. i' P/ x; x8 G$ g0 V' d" `3 W
 for (i=k; i<=n-1; i++)+ r7 i" E/ y1 G" N" c
 { v=(k-1)*n+i;
( h# b0 ?& S' W, m3 w d=d-a[v]*d[i*m+j];
( E, R( X& S& Q2 m7 m H/ U d }8 h# i3 r5 u" U) Q$ g
 v=(k-1)*n+k-1;
+ O. g7 a" |8 g5 E0 ~2 a* a d=d/a[v];
. g+ q3 Y. N! u1 P- X1 D3 \) x }
/ O- @" \( H _5 ? }
& a G0 `( F P return(2);4 J3 g* c4 |' B" k, `* J, d% a5 ]
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)
 {1 j$ G5 u! ~4 y6 \; x0 v int n;. Q; a% | w$ p' T; J
 double x0,h,t,y[];
! O: p0 e& w4 N& D C { int i,j,k,m;
2 U+ Z2 G: F6 d- B double z,s,xi,xj;
8 a& w' B& u6 E: Z float p,q;& e: {, D1 z( K& Z% b
z=0.0;% m* b' b o+ T6 i- H
if (n<1) return(z);/ w' v9 s9 X6 O0 O. w$ y& d
if (n==1) { z=y[0]; return(z);}2 o5 w# A+ |. \$ Y$ o& h. B0 @
if (n==2)! J. @- G; j* r
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;8 @9 l% s/ M S( o! U0 A$ }
 return(z);3 X8 C% ?% W" Y ^& Y" E( N: G& o
 }$ Z3 o0 h! ^3 r+ k& o) S
if (t>x0)
$ g9 x0 q" S3 ?6 q# Q0 x9 q. F5 L' ~9 @ { p=(t-x0)/h; i=(int)p; q=(float)i;* m; d, Y( R0 K, k4 v6 f
 if (p>q) i=i+1;/ m1 E9 s9 q; K3 t
 }
! N. H3 _9 k& Z else i=0;
0 n4 F! h" |7 { g1 j9 I( X k=i-4;
& H. s7 R2 f* L1 h( T% f if (k<0) k=0;
9 J+ R( ?0 k( y- K+ T* S% V m=i+3;
/ \9 r% U; @4 [3 m if (m>n-1) m=n-1;
3 s: z5 e* x" F K! H- e7 ]3 m& v for (i=k;i<=m;i++)4 @/ Q6 B) g* i
 { s=1.0; xi=x0+i*h;, [+ \4 K: m& f1 Q8 x+ b/ x! Q6 c/ B
 for (j=k; j<=m; j++) W- ^& V/ b/ b# w& R
 if (j!=i)
1 C- G7 ]+ k- s( e4 \/ H% j { xj=x0+j*h;
 {# q: u% R' ~8 J9 }# y s=s*(t-xj)/(xi-xj);0 i: Y* t3 U! y
 }' u1 d( Z$ X8 a1 k) j
 z=z+s*y;
7 H2 z) D) A& ^) L' v }1 L# M R+ t8 l
return(z);4 N2 c/ P' w: \ M
 }; Z3 l- ?2 q0 K1 H
向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"
/ x8 {# |+ @0 V4 J% t% `; T8 y! W void hpir1(x,y,n,a,m,dt)
7 J$ D8 P9 v# u int n,m;
" x7 g) Q6 l; S) f) H0 } double x[],y[],a[],dt[];
. B/ @5 l- Y, F( d { int i,j,k;
$ Y6 R0 o- F5 S& z& D double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
! [: A$ o& A. X' a/ x7 ~9 l for (i=0; i<=m-1; i++) a=0.0;
4 ]! R! ^( l2 }! ]1 Z if (m>n) m=n;0 o% U3 e% w8 p
if (m>20) m=20;/ s7 j/ L% M: j9 q4 @+ d! L
z=0.0;" n h6 y" ~- z4 h$ B9 k0 a
for (i=0; i<=n-1; i++) z=z+x/(1.0*n);9 V1 l0 O7 z; s- l
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;0 |6 k" W% z$ @9 `: _( v
for (i=0; i<=n-1; i++)
/ D' n9 |. m: B9 _ { p=p+(x-z); c=c+y;}
. K$ ~% n) f/ M/ e* T& l: ] c=c/d1; p=p/d1;' U. ~& K, i% y( X
a[0]=c*b[0];
4 c$ u& ^( _6 Z* D/ A if (m>1)
1 T3 ~* w/ P& E+ }1 r3 g { t[1]=1.0; t[0]=-p;9 U& A7 L4 ~1 r" }# y7 q
 d2=0.0; c=0.0; g=0.0;& |/ H- V _4 Z* y/ l; _; O
 for (i=0; i<=n-1; i++)/ F9 q' n4 a2 U- ?; T o
 { q=x-z-p; d2=d2+q*q;& P3 i) ?. N, _9 z% l" S$ ]
 c=c+y*q;
+ Y+ r0 D; t' J4 f( u g=g+(x-z)*q*q;
) E. y0 `3 G& R4 g7 I f! @ }
: n, z: R" ~* C- x! A% { p+ [ c=c/d2; p=g/d2; q=d2/d1;5 J. u2 L+ @/ Q- [7 M1 [& @
 d1=d2;8 R" I9 K3 N' |/ W5 L6 \3 x5 y
 a[1]=c*t[1]; a[0]=c*t[0]+a[0];
7 ? c2 k9 I* C; J2 G3 U }, [ t4 v# C7 K6 g {) l$ P
for (j=2; j<=m-1; j++)
+ w: u8 q; K& }0 P' n { s[j]=t[j-1];/ c1 g3 n3 t3 [6 D; f3 r) ?
 s[j-1]=-p*t[j-1]+t[j-2];5 Q2 E1 R! H- A( O5 z8 J9 @
 if (j>=3): n+ }% I P. Q9 O) F: G$ a+ C- ~
 for (k=j-2; k>=1; k--)1 w: M, F0 d" w2 U
 s[k]=-p*t[k]+t[k-1]-q*b[k];1 g0 [8 {; A1 Z2 J0 ~' x" m7 D
 s[0]=-p*t[0]-q*b[0];2 i N) m* Z4 \/ W4 o+ C
 d2=0.0; c=0.0; g=0.0;" i' l1 x; |7 W. R* j/ t2 h
 for (i=0; i<=n-1; i++)
9 u% `+ T& w- E { q=s[j];0 t. j/ ~3 w* y
 for (k=j-1; k>=0; k--)
! A: A+ X3 F8 }3 i5 z3 d y q=q*(x-z)+s[k];
5 Z$ g2 F' C# j* [" b d2=d2+q*q; c=c+y*q;9 I3 j5 R: z, H- B8 j( \* q% ~
 g=g+(x-z)*q*q;
5 n1 b8 ?% U! v5 t. h, b0 @( s+ q }
& |( [/ ^- V' I8 K c=c/d2; p=g/d2; q=d2/d1;* B3 U2 L% v& p p+ }- ~6 {
 d1=d2;* J- L' H S& J! y) D4 o# C' \
 a[j]=c*s[j]; t[j]=s[j];* t. M- v, ?7 B* f4 Z
 for (k=j-1; k>=0; k--)1 h0 T* b& o5 Y% @
 { a[k]=c*s[k]+a[k];6 X# s; _* B n9 |# m F% H9 Q. h
 b[k]=t[k]; t[k]=s[k];
6 [9 r9 E, o& I) r+ E2 X, n }& M- G& ]; w; n9 q0 {; A
 }, G% H: o9 X" S. V+ ] e8 b
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;* k5 a& d$ j$ i. U" u0 H6 T
for (i=0; i<=n-1; i++)4 w! q. i6 Y* z) J
 { q=a[m-1];# h- v+ e% i Y7 `5 @6 u: k/ M# u
 for (k=m-2; k>=0; k--)0 @/ M5 F" m8 \' f+ h& |+ p, Z
 q=a[k]+q*(x-z);7 r+ e0 k6 R7 ?' p& L
 p=q-y;
, S/ G$ u# C9 m1 D if (fabs(p)>dt[2]) dt[2]=fabs(p);$ ?& b6 m) Y4 @% S+ ^
 dt[0]=dt[0]+p*p;0 l6 c5 O$ b3 G4 l/ ]) N* u
 dt[1]=dt[1]+fabs(p);# [+ j4 `: S) q6 b4 G9 J
 }1 A2 P) e# j; p# @
return;
5 h4 K+ ]$ I* i& P, ] |6 @ }<>龙贝格积分法<>#include "math.h"
# t V# C) `3 H' c* ?/ h double fromb(a,b,eps)2 R2 ? b; x1 j+ \
 double a,b,eps;
3 A9 K ^. C R5 K { extern double frombf();* ?- Z. h1 I" b
int m,n,i,k;1 ^' f q/ q, f+ H3 v% n4 R
double y[10],h,ep,p,x,s,q;
7 [6 J2 t! C5 s) U! G. d1 w h=b-a;
8 E# R+ y2 i" ?) I, p y[0]=h*(frombf(a)+frombf(b))/2.0;
( h! T4 N1 b; o! e7 Q* z& x m=1; n=1; ep=eps+1.0;
3 }+ n$ c q. h9 g Z5 M% P. k while ((ep>=eps)&&(m<=9)). s) z) Y6 }( C2 \8 r
 { p=0.0;
7 k9 p; ]' @+ q for (i=0;i<=n-1;i++)
8 K( H. i+ q& { { x=a+(i+0.5)*h;$ f- f$ X/ ~/ r" D
 p=p+frombf(x);
, P# m5 e3 i" r7 I, k8 ?3 [* T% E }
/ B; {. }) A! ~$ z' [. l+ K9 H p=(y[0]+h*p)/2.0;
% P- f) k- z l9 h( Q$ K' h1 N s=1.0;
3 M% ?! P9 o9 Y* W for (k=1;k<=m;k++)9 t9 j1 w( S- w2 v4 M9 M K( v. c
 { s=4.0*s;* x& W9 _, @( r4 B, z
 q=(s*p-y[k-1])/(s-1.0);0 _7 H% k( Z6 K- G3 l
 y[k-1]=p; p=q;
! s6 U, _" u9 s- P, h }% J/ k. L8 z/ d S) ^- Q
 ep=fabs(q-y[m-1]);
$ q: z8 f$ F* Z, N7 T m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
 K% `8 S: k1 W }
8 b9 M# k/ L6 B' G. A5 B3 j3 B return(q);1 `& n$ M. g: T7 v" E
 }<>呵呵希望对你有用！！