查看: 5958|回复: 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;
, W. m+ u9 H9 H* N" \& r3 O }
$ ^# m7 W, I+ s }
2 O/ y: O5 ? r7 U! J- S u=(n-1)*n+n-1;* l$ W; f- d! ]" k) b9 m
d=ar*ar+ai*ai;
7 O" P8 Q* @' T* B: k if (d+1.0==1.0)
0 C r6 }: o3 y { free(js); printf("err**fail\n");
& K2 X1 `! M& m5 b return(0);2 K+ [, V# r4 Q6 u. n, M! e# Y) M V
 }3 G) {: R( P& y$ ]( i# h/ H
p=ar*br[n-1]; q=-ai*bi[n-1];
2 r( N; b7 v) |, I1 N+ D& D: P s=(ar-ai)*(br[n-1]+bi[n-1]);) W- u+ H8 y3 o. B. f7 |5 u9 G+ ?
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;) `5 A! C/ Z0 R" @" L L' G
for (i=n-2;i>=0;i--)
' c/ B$ _' H, `* q4 [# u for (j=i+1;j<=n-1;j++)' K) H- e7 Q: i5 \+ P2 F6 Z9 Q
 { u=i*n+j;8 M. q5 W( L2 {. g
 p=ar*br[j]; q=ai*bi[j];
+ D! A* V1 g: b1 \8 O s=(ar+ai)*(br[j]+bi[j]);" `/ m; q6 N' S! p
 br=br-p+q;; e" Q8 M! \: t0 y" \& J3 n+ K
 bi=bi-s+p+q;
1 b/ |5 f" l; r. |% @4 M) D! e }0 ~& V% u; ^+ C/ c: z6 x% H
js[n-1]=n-1;
; T- }3 y" k4 J# n1 x2 ^2 S for (k=n-1;k>=0;k--)
7 Z6 f% |" I4 N: b+ y4 T0 i) P if (js[k]!=k)
+ Y4 p: u: P, c" Q$ } { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
1 d' y/ l- [% X% r) @. j p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
+ g9 S$ c: Y/ H* Q& U8 Q }! l8 ]9 ]; K% b1 E8 h& s! t% r
free(js);: x, w% E, M+ e3 R/ l
return(1);* C' K! h0 C( E, }- k; x
 }<>平方根法<>#include "math.h"
9 |3 W& h4 V. H" A8 J" h #include "stdio.h"
 A# S' P, C. r' j- i. r int achol(a,n,m,d)$ [3 X) O9 Q% F) F
 int n,m;3 |. m2 F" N% X0 u( v
 double a[],d[];8 ~7 m$ S/ j/ S% _: \! c; x. O
 { int i,j,k,u,v;
" _" Z4 M7 h5 v! {' s6 ~! V if ((a[0]+1.0==1.0)||(a[0]<0.0))$ \3 z: Y* h) G' v. x" ?
 { printf("fail\n"); return(-2);}" P p/ U# Z. Q5 J8 l
a[0]=sqrt(a[0]);
- U7 I. I6 N* j& P! z9 o. H for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];5 w; X! [) N5 R- M- g& N
for (i=1; i<=n-1; i++)
( V4 A: w+ Q+ E p { u=i*n+i;
$ c7 j) |, I# p x for (j=1; j<=i; j++)# ~$ j6 K. D+ c& K
 { v=(j-1)*n+i;" _) C4 K6 n& v7 s' M: t) H3 B
 a=a-a[v]*a[v];
2 Y9 o c7 I; \- j }
9 `9 U) R) D3 V! r2 C, @9 l' N if ((a+1.0==1.0)||(a<0.0))! w, j" h$ X9 H# v
 { printf("fail\n"); return(-2);}; z; k* w* ]6 \8 I% ^ J
 a=sqrt(a);
* F; r5 q& ^: O, p% N if (i!=(n-1))6 {0 \/ [3 d) j+ E3 `
 { for (j=i+1; j<=n-1; j++)
) T7 t& H: q# ?8 r# S2 t: ^9 b$ i) _ { v=i*n+j;( o4 @9 I! u6 V8 v3 f, O1 c! V6 g
 for (k=1; k<=i; k++)
0 Z) y& L& T$ ^( @9 y. o) b# M. ~ a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];/ h9 o- l( C, ]( \( b9 o6 |7 V0 p
 a[v]=a[v]/a;7 m, Q6 K: T0 S- @9 q$ `
 }. u; I* c _" T4 {/ S
 }5 \1 t# u# j) b( m
 } C8 T8 O3 F) W! E" \' a% P8 L
for (j=0; j<=m-1; j++)3 k; D$ l. E$ i r' {% t
 { d[j]=d[j]/a[0];$ p5 S+ O" J6 K. H
 for (i=1; i<=n-1; i++)5 E$ }0 k8 t8 _, N* z& ?
 { u=i*n+i; v=i*m+j;
% D& I, Y7 X, Z! Q6 x) o for (k=1; k<=i; k++)0 O; l) m/ y3 f
 d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];% q8 b8 c9 Q0 [1 K- G" B
 d[v]=d[v]/a;
5 C% {" x" g) z. R' n }6 o4 A8 P0 U' j- k+ L
 }
3 v6 r$ W/ N- ^( G5 s for (j=0; j<=m-1; j++)
# I2 ?; z( [9 a6 S0 c" i! G2 q; w { u=(n-1)*m+j;
7 e F- |' g7 R V$ c d=d/a[n*n-1];. p: ^# W( P4 @: I* ~
 for (k=n-1; k>=1; k--)
: I1 p% N X: j& E { u=(k-1)*m+j;8 C% k3 B P6 l
 for (i=k; i<=n-1; i++)( Z- \8 ]! J. v% S
 { v=(k-1)*n+i;
- D2 ?- l. ?& v) L& O; x0 { d=d-a[v]*d[i*m+j];! r/ g9 y+ T0 G7 A2 Q ~! o/ v0 O/ R
 }
; w; d. }+ s( q. r& |: S v=(k-1)*n+k-1;
6 T% s+ \ F0 X3 y* E) d; E d=d/a[v];
& _7 C/ M: P) h& J# W }. n; {* u9 v B n! r
 }
1 g4 M2 `5 Y! j& \ return(2);( d* }$ t% l6 M
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)
9 k8 Z" |! o8 I# w3 a Z% V int n;, J5 P/ ?" H" _) t/ U1 u: D' H% x$ b
 double x0,h,t,y[];2 ^4 M" @$ u- H& V7 l5 C9 t
 { int i,j,k,m;
8 D- j, |- }( s5 Q- w double z,s,xi,xj;% i% \$ \8 `# ^+ E6 p4 d4 ^" c
float p,q;! i( d) F+ h' l9 |, T& ]
z=0.0;
/ Q# Z" Z& V8 m' X! v- v if (n<1) return(z);
+ P+ ~6 @3 x2 o if (n==1) { z=y[0]; return(z);}
& r/ V; b8 ~* w1 t7 |$ M4 [ if (n==2)7 f2 X$ w5 R5 ^0 o' l
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;6 M* L) k I1 K( \: a$ o7 M
 return(z);
/ m+ T9 m; A! O2 ` }
# t0 P3 E9 _6 K: L+ }- N if (t>x0)
: o8 |6 V3 c$ ]1 O6 [( e( l% ` { p=(t-x0)/h; i=(int)p; q=(float)i;( K# D- N% A3 q' k$ o
 if (p>q) i=i+1;4 _+ h- ], f9 v6 b/ y$ Q
 }' t# ^3 `, L: E+ W
else i=0;
; v1 G6 v( j Z9 O X. e- E: d( C k=i-4;. b) x% E9 B, R
if (k<0) k=0;
2 o, q \3 |( p3 M( @ m=i+3;5 d* ]6 }# K1 c
if (m>n-1) m=n-1;" v" T- q+ T: C' D% Y
for (i=k;i<=m;i++)& T# L+ K+ R, N6 X- B" s" C8 k
 { s=1.0; xi=x0+i*h;4 p$ q/ M; z6 d: h/ N+ m
 for (j=k; j<=m; j++)6 v8 o, U. ]. E r
 if (j!=i)1 W" `3 Z c4 k2 z! q6 n9 K
 { xj=x0+j*h;
1 l- Z7 U( U, ]9 i s=s*(t-xj)/(xi-xj);0 g/ I, o: L) |9 V
 }$ z2 {9 s+ F2 b
 z=z+s*y;
' H6 F6 H1 C3 `/ h5 p9 @: o* Z; T: r# K }; y' z0 W& D2 V1 t% _! ]! J
return(z);. S% |! k0 }7 J0 p8 k
 }$ h# K; X0 |% H5 r' B+ E
向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h") g! a9 v! `6 t" F
 void hpir1(x,y,n,a,m,dt)
2 ~- ~ q: P! K) Y9 { int n,m;
; @- C& W Q' x* t9 S# L double x[],y[],a[],dt[];
9 A2 m( W* z5 W1 J1 { { int i,j,k;
3 }4 L# \* M0 ^0 K; ^8 p5 A7 P double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
" A) f0 c# ?% Q/ C o. |# R for (i=0; i<=m-1; i++) a=0.0;
' \6 z- X* K& Q; M: v+ d' f" \ if (m>n) m=n;
' N$ W7 p* }$ B8 N& E* k8 O( s" a if (m>20) m=20;
6 A" G4 a8 ?7 d( a- ` z=0.0;
: G% C1 p$ x. { for (i=0; i<=n-1; i++) z=z+x/(1.0*n);
 r: N" u5 R1 l& _5 S1 N b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
1 Y! F, D% U& t2 T4 h) E" x) r for (i=0; i<=n-1; i++)0 @2 x8 T# i) {: |
 { p=p+(x-z); c=c+y;}
) X; r- i6 `1 y$ H! E! c, c) U c=c/d1; p=p/d1;2 F! d2 a3 K; k! D0 `
a[0]=c*b[0];* S; w4 x" j' P8 t2 G
if (m>1)
4 K* y, w. V( H5 w! Z { t[1]=1.0; t[0]=-p;. M. H! L+ ~7 A+ y) M3 ^0 C
 d2=0.0; c=0.0; g=0.0;
# p# y4 `. v1 M1 m2 b R for (i=0; i<=n-1; i++)% [$ @4 g" L/ b/ ^) b/ y9 B$ n3 M
 { q=x-z-p; d2=d2+q*q;
* L) I" f5 Y0 D! m c=c+y*q;5 |. W2 m X8 Z( m; H' j3 _; q+ I
 g=g+(x-z)*q*q;7 b4 e& ]3 F. |6 c6 v1 w
 }( d2 O2 P5 B( w+ ?# v
 c=c/d2; p=g/d2; q=d2/d1;! f3 n, m9 i5 z$ ?
 d1=d2;$ Z: R/ d" K. T
 a[1]=c*t[1]; a[0]=c*t[0]+a[0];* d, f* j* g7 _5 |
 }
 B6 }# i" y0 y$ ~ for (j=2; j<=m-1; j++)
 U$ T2 g; X q( \; g' | K { s[j]=t[j-1];. ^" [% w5 l1 Z' O( M
 s[j-1]=-p*t[j-1]+t[j-2];
/ B3 ~5 l& ]: _2 [. U if (j>=3)0 r' r H' U! ~3 S
 for (k=j-2; k>=1; k--)! i; O. L; y' W2 V
 s[k]=-p*t[k]+t[k-1]-q*b[k];; \' {* l. j+ P8 c
 s[0]=-p*t[0]-q*b[0];
; u. G* h! \" e* Z! Y( J. u6 ^7 x d2=0.0; c=0.0; g=0.0;5 e+ u& Z p' s
 for (i=0; i<=n-1; i++)
' `, V( M' a, V2 m& x { q=s[j];
+ h" W! v% E1 F& q for (k=j-1; k>=0; k--)) e+ S# n7 o8 F' J6 G
 q=q*(x-z)+s[k];
( I* C# M* A+ c6 T8 T) t2 G" G d2=d2+q*q; c=c+y*q;! I8 M% J% ?+ U8 g) ]3 j m2 N
 g=g+(x-z)*q*q;% r5 K1 x* i9 L! X& S
 }' t8 O4 Q) c* \/ ]5 z
 c=c/d2; p=g/d2; q=d2/d1;: a1 m7 I' C# ]+ y* {4 s; O9 O
 d1=d2;* q* n: E9 M' A2 M, z; E
 a[j]=c*s[j]; t[j]=s[j];* j/ X9 [- H$ f0 W7 ]
 for (k=j-1; k>=0; k--)3 C8 H4 U' T, S ]. f3 ^' J4 t
 { a[k]=c*s[k]+a[k];/ H: \- J. Z0 z
 b[k]=t[k]; t[k]=s[k];+ P# E5 {' F% e* ~8 Q' `7 j
 }2 g2 O4 J5 {1 N& f6 s
 }: m+ j( W- K* I+ A1 w( W7 B
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
5 ?1 f; S, }; f+ U- Z8 q- { for (i=0; i<=n-1; i++)
/ G p: Z5 ?0 l1 l { q=a[m-1];
, C; i% [( F$ Y4 ^$ H& J9 z- ~ for (k=m-2; k>=0; k--)
6 G; X# x2 f) a8 x) G- H q=a[k]+q*(x-z);9 } r1 x0 q. O! P4 j
 p=q-y;$ e; Z* @( o, R/ o/ {
 if (fabs(p)>dt[2]) dt[2]=fabs(p);
; b) F1 \8 K9 \$ x5 I dt[0]=dt[0]+p*p;
4 q- B& K# N* ]: B/ h" \7 ~( C$ Y dt[1]=dt[1]+fabs(p);, a! `8 K; y# h( F
 }6 Z+ S, R% o: h6 M& b* X
return;
. A d5 i3 n% Y( Z* s7 T }<>龙贝格积分法<>#include "math.h"
: e& b. [# m7 H double fromb(a,b,eps)
( d& e+ m0 D, Q. W* |7 R& A7 n double a,b,eps;
; y. x' |1 z7 v: W: G { extern double frombf();
4 z7 A5 H0 X8 |6 |! o g- s, _: M int m,n,i,k;
5 P, l: i5 S+ { double y[10],h,ep,p,x,s,q;( \* s* }5 Z! K4 t1 q
h=b-a;2 q! T+ \, ?3 Z$ b
y[0]=h*(frombf(a)+frombf(b))/2.0;' F m' b9 h4 f5 k* b
m=1; n=1; ep=eps+1.0;) W2 C' e! |9 U' Q; w+ M5 p
while ((ep>=eps)&&(m<=9))1 r+ l2 K3 j* s0 A% I$ P: z5 U0 V
 { p=0.0;
/ M/ f0 U' W' `2 R: z% c for (i=0;i<=n-1;i++)
1 E: @! b) W+ n9 H7 T { x=a+(i+0.5)*h;
5 g; h) |6 c2 W% \ p=p+frombf(x);
# u% m1 O7 @" c( P# B! z4 { }. g- |- k* Q1 H3 y, D: ?& [
 p=(y[0]+h*p)/2.0;
0 ?4 f7 |/ i7 J7 h' o0 j3 g s=1.0;& q( A0 `3 O0 v7 Y, ]7 X/ A4 o
 for (k=1;k<=m;k++)/ U- E. Y- ]* H9 Q7 d4 P
 { s=4.0*s;* r$ j7 s; H. }$ [
 q=(s*p-y[k-1])/(s-1.0);/ L% ^4 ?/ |8 f2 d/ t
 y[k-1]=p; p=q;7 K% ` Z0 g2 a; M6 f9 Y
 }
* o! s+ a5 R x; y& X3 F6 v ep=fabs(q-y[m-1]);
4 S1 p8 q: J7 w1 A) W* @ q/ B m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
, l, C& N# e3 x }
2 |6 V$ c4 m" R return(q);, j) H* y6 i0 a C, M
 }<>呵呵希望对你有用！！