查看: 5963|回复: 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;; C9 l5 `3 {2 x8 E* R1 p. K
 }
+ W9 e' d0 ~+ C7 @3 q4 u! u }# Z% T* A3 Z# g5 B
u=(n-1)*n+n-1;
7 o, g; i0 W2 C d=ar*ar+ai*ai;' P# \ f$ o, r X
if (d+1.0==1.0)0 x; V! H! ^- U( G
 { free(js); printf("err**fail\n");2 Y* I% {- J6 T: M# a* u( I
 return(0);! _% l; M+ ? T9 X$ d& x
 }
/ a7 h( O3 @7 {2 Y' d p=ar*br[n-1]; q=-ai*bi[n-1];4 O4 r3 o4 M9 F/ h! ?
s=(ar-ai)*(br[n-1]+bi[n-1]);
0 t5 i) [* n" u5 B br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
4 Q& P5 T3 y5 A$ j$ w for (i=n-2;i>=0;i--)
/ ~- X/ j; s- ] C# r. }4 ]* h for (j=i+1;j<=n-1;j++)
0 e- ^5 j& ]# J' m. t { u=i*n+j;! u0 N9 A3 n% {5 X- [1 Y
 p=ar*br[j]; q=ai*bi[j];
4 o* a) r% q# ?6 b, I) S8 A s=(ar+ai)*(br[j]+bi[j]);
/ Y- S: a# l& h# b br=br-p+q;
% Y2 A: ^4 h3 e9 ^ bi=bi-s+p+q;0 f4 {* j ]/ U' j
 }
/ n6 J# g' p/ h; K9 R js[n-1]=n-1;& V' n! g) ^( K" ^+ @
for (k=n-1;k>=0;k--)
7 E& u9 b& D; v' h7 C) | if (js[k]!=k)
7 A+ K& ]2 G) M { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
. _9 O# ?1 z! m5 V- E p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
) r# L0 s4 }% P0 b1 I0 f1 w; s! ` }
+ j! q- r0 l/ x3 b9 c- ?/ h; U" Y free(js);% n/ m" }5 Y$ R& _: ~
return(1);
. S; h3 |( T5 r9 O: H }<>平方根法<>#include "math.h"
* [5 V2 h# O5 x+ P, e #include "stdio.h"1 W# H. r |1 ]* }! z- Z' |" Z
 int achol(a,n,m,d)7 \: N* k( c' [( g
 int n,m;1 H( k( v3 E. T
 double a[],d[]; k B5 ?3 v1 L% r4 W; M
 { int i,j,k,u,v;
) H& U2 N2 j! \. z% r, N! x v if ((a[0]+1.0==1.0)||(a[0]<0.0))# F; ~# q8 \' `+ ^6 }3 j
 { printf("fail\n"); return(-2);}
; N2 a5 P# V" y a[0]=sqrt(a[0]);1 ?, R4 P c- i4 c$ X$ r
for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
* B+ [ F2 N. d' K8 p# a for (i=1; i<=n-1; i++)
2 J7 Y% _: _. S4 q( w6 V { u=i*n+i;0 H S& _; Z, `, O' E0 `2 s7 @
 for (j=1; j<=i; j++)# Z) q& Y/ A9 G( K6 S
 { v=(j-1)*n+i;' z6 i* y2 W# S$ w, r
 a=a-a[v]*a[v];
 b9 J" P9 ?2 m: i1 @. m( x }
 a: x4 G7 h. }4 p; M if ((a+1.0==1.0)||(a<0.0))8 T Y, p, j. f" Y% y' g- k7 _
 { printf("fail\n"); return(-2);}, ]) V. V. a6 M/ a* D; G
 a=sqrt(a);% ^. {0 t# I! [! H
 if (i!=(n-1))
$ V7 q0 M/ T0 s: ^* ?4 } { for (j=i+1; j<=n-1; j++)' Q! ]9 {8 A: p5 a$ G
 { v=i*n+j;/ u" U% H2 O7 j* s
 for (k=1; k<=i; k++)0 h$ M0 z$ K8 `. t
 a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];: \5 i; H b* `8 w0 u$ V8 ?1 R
 a[v]=a[v]/a;1 h# b1 S" ~+ q) ~+ _
 }
- @4 o3 q% B V }" a+ I: x% F# {& t5 H2 M0 G6 L. J
 }
& Y2 c& @ `( g# h for (j=0; j<=m-1; j++). L% A/ a1 u4 c3 E3 e" r2 M4 M% z
 { d[j]=d[j]/a[0];) t( n* T( f0 X" s7 i: [
 for (i=1; i<=n-1; i++)& E7 Z M/ h9 b% q; T ~. U
 { u=i*n+i; v=i*m+j;
: |% w l9 g {- p& P6 A for (k=1; k<=i; k++)
7 k& p/ d5 }, }+ }6 j d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
6 } L! h) v3 E" m- N- d& R d[v]=d[v]/a;
5 R3 t/ o' Z- X r }$ X; T% w4 }9 A2 Q K; r7 Q
 }5 {4 I% k5 M3 K. Z& G$ J- T6 s
for (j=0; j<=m-1; j++)
3 h! e: P! ?6 } { u=(n-1)*m+j;% s( C5 h& T! |4 a
 d=d/a[n*n-1];: s3 f2 _- N6 C$ w# L( k
 for (k=n-1; k>=1; k--)/ @% F6 ^8 a+ ]8 N
 { u=(k-1)*m+j;7 A, @! o" Z& U3 ~' }
 for (i=k; i<=n-1; i++)
, @5 m* X" ]+ D! d+ n { v=(k-1)*n+i;
6 b( |3 {9 \. y$ ^8 i o2 T d=d-a[v]*d[i*m+j];
5 o% M5 @1 u( B }" g, {3 K, r3 A$ m% s
 v=(k-1)*n+k-1;+ G3 X4 X( X, [, N! ]9 |: O4 r" ^
 d=d/a[v];2 ^" h- P' A8 F2 S6 ^1 x
 }
$ g! a2 M1 y9 q+ r }
/ C( V+ L& n3 b' r4 D- O5 c return(2);4 c2 U* Y4 z/ @# c
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)2 a. Z( I- k, X( H8 q$ \% P
 int n;
4 }! J: y6 O8 C& g& t& _2 p9 c- w$ F+ [ double x0,h,t,y[];
5 W; k9 F) ]: s8 { { int i,j,k,m;
8 f- \. x U3 J' K& C$ n+ m double z,s,xi,xj;/ q( a: C: L7 P
float p,q;3 ~" {( c# }3 S( P
z=0.0;
1 z; u1 v3 W8 K S5 n if (n<1) return(z);2 k/ H% _" a7 a3 _+ u+ U9 @
if (n==1) { z=y[0]; return(z);}
& v* R* S) r! M( \; Q( I if (n==2)
* r. ^0 J2 A3 W; M w$ e# h { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;' o5 v# e7 p+ [: R& j3 |
 return(z);# [7 C( x& B" {, q; e
 }
; ~% T! ~: E+ D+ A! i* ` p. M if (t>x0)( z, ]* X0 q& u1 v
 { p=(t-x0)/h; i=(int)p; q=(float)i;
$ U4 e+ |: R$ h7 D) e+ X if (p>q) i=i+1;! D. i# t) i5 x* \7 F
 }
! d1 H( ]; b! p* R7 ^2 r2 s else i=0;
3 R1 i' B4 e' q! s* K4 E! Q! V k=i-4;
7 ^% D. g! Y" [ if (k<0) k=0;% r8 f6 w6 K% r6 _
m=i+3;& g; Y* c0 ^4 ?2 K+ r/ K: B
if (m>n-1) m=n-1; h" {1 G3 K. \8 W
for (i=k;i<=m;i++)1 V, M$ J- a& l8 g
 { s=1.0; xi=x0+i*h;
* E. a: V0 u$ u/ G6 Z; s9 G7 c for (j=k; j<=m; j++)
 k' x6 L! n" }% p: `1 H: { if (j!=i), z* }8 y+ H4 O) }
 { xj=x0+j*h;
& t3 O6 B3 I5 G; U s=s*(t-xj)/(xi-xj);" ] }6 l+ [9 V# |* y1 V% _, U4 G
 }, Z: u6 ]- K+ T9 D2 }! g2 G6 v
 z=z+s*y;
$ l0 \9 n8 z' K }
$ x/ `+ ~6 u) j0 j0 ]# j$ n$ s3 m return(z);0 y0 j9 [3 x! e
 }" c2 V5 @% R% P3 ]* L
向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"
, _7 W% p% M" @; y7 g. l void hpir1(x,y,n,a,m,dt)
1 c9 Q7 P5 Z. H; Y int n,m; @+ g" o% G" ~" Y5 B
 double x[],y[],a[],dt[];: s0 X1 W7 v6 ?" L) b' v! X8 f
 { int i,j,k;
! W6 V% B" }. e double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
2 i3 e( a. o, B# Q/ X3 C for (i=0; i<=m-1; i++) a=0.0;
9 u; d$ W! M$ V! c3 ?0 |* [8 { if (m>n) m=n;* W0 V5 [9 M+ n0 ~) Z7 O6 q4 f
if (m>20) m=20;$ {4 k/ o, y5 f E' i
z=0.0;/ k, i! r& E9 O
for (i=0; i<=n-1; i++) z=z+x/(1.0*n);
1 X7 A* N! {, a& d6 e, f b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;+ G9 I0 ~8 M! }5 k
for (i=0; i<=n-1; i++)
$ T+ J. @( ]$ T$ t { p=p+(x-z); c=c+y;}
! C: z4 W+ _1 B8 ?* L: k- v0 v& A) G c=c/d1; p=p/d1;6 }) z: _ L5 b. X
a[0]=c*b[0];
3 C1 G" t2 [8 p7 j' H& j1 Q if (m>1)$ z O0 c1 y1 Q8 J
 { t[1]=1.0; t[0]=-p;( m" D7 f( o9 N7 J+ l+ G6 `. ~6 D
 d2=0.0; c=0.0; g=0.0;/ X2 H; a X3 ^
 for (i=0; i<=n-1; i++) y5 `$ i' F8 D: w/ B- \
 { q=x-z-p; d2=d2+q*q;6 u1 r9 T* ~# z. l" Z+ X
 c=c+y*q;/ l1 p1 A0 S3 c+ X$ `
 g=g+(x-z)*q*q;2 m6 n6 T' D+ V# m$ z, F* J
 }
5 M o3 l) h9 F/ O+ J1 p* N c=c/d2; p=g/d2; q=d2/d1;
: p' S/ s: E# S8 t d1=d2;
' g8 V0 W4 Y8 i a[1]=c*t[1]; a[0]=c*t[0]+a[0];1 `$ k( v8 _( o: J0 K5 r$ L
 }
8 a: t8 W" S/ _: p" z- @& \ for (j=2; j<=m-1; j++)
! z9 i. Z) B* A { s[j]=t[j-1];- w9 {3 \2 m" n
 s[j-1]=-p*t[j-1]+t[j-2];) s6 H. L! h& R# K! G# a" s( j0 v
 if (j>=3)$ b. f$ F4 w% u9 G: m; H/ c+ }3 e
 for (k=j-2; k>=1; k--)' T; q0 K, q% F
 s[k]=-p*t[k]+t[k-1]-q*b[k];8 h/ G/ ^% |% g$ N9 p
 s[0]=-p*t[0]-q*b[0];6 B) r& n9 }/ `, H' U" B5 e
 d2=0.0; c=0.0; g=0.0;9 Y( a. L5 O; J
 for (i=0; i<=n-1; i++). D! e8 h8 M, K+ o5 O
 { q=s[j];
6 H8 U" T7 x& X3 O4 k8 C for (k=j-1; k>=0; k--)
7 t# Z' u% N3 s; G q=q*(x-z)+s[k];
$ A& r* c0 H" A w5 ? d2=d2+q*q; c=c+y*q;
' ?7 @/ W# X6 b g=g+(x-z)*q*q;
; W8 p+ e8 I# r% I/ E% v }' A! T. \2 g% \" h* U& U* S
 c=c/d2; p=g/d2; q=d2/d1;8 z* F4 g1 @2 A/ |& V
 d1=d2;
4 Q5 q% g- y$ A# z/ }$ H$ [5 u a[j]=c*s[j]; t[j]=s[j];
' A& E @ K2 ?, K. w3 {# W: ^ | for (k=j-1; k>=0; k--)$ k9 C# R }6 v4 \
 { a[k]=c*s[k]+a[k];! ` f6 s# `/ D6 p
 b[k]=t[k]; t[k]=s[k];
- A, N2 J' Z+ ~2 ]0 n6 U0 b }( u. r- x% E( e: z
 }4 E; Y5 u+ d2 g
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;0 |( g, L* I: L3 l g! ?6 j4 u
for (i=0; i<=n-1; i++)
: S4 C, E) a2 t { q=a[m-1];
( o: e+ `& q) x3 {$ g' M& ]* T) |9 P for (k=m-2; k>=0; k--)
+ a0 X6 d4 Z1 N. M$ ~ q=a[k]+q*(x-z);
6 J* f/ F* `, ~2 s# W p=q-y;, }; @* h# c% q/ V7 o5 B4 Q6 U
 if (fabs(p)>dt[2]) dt[2]=fabs(p);1 ]) l+ N; ^1 ^% f; N' P
 dt[0]=dt[0]+p*p;/ y4 V' L5 ^* O- y- ~( l
 dt[1]=dt[1]+fabs(p);
! C% R& D3 R4 C* x4 T }) K9 Y. G9 w5 E* F' R" j8 ^
return;) j; ^- X; D/ t0 |8 S! L2 }
 }<>龙贝格积分法<>#include "math.h"' K1 M$ ?9 X% c% U& v5 @& ]6 L& \6 F
 double fromb(a,b,eps)
0 t2 M I5 \' w" Y5 ]) K! \* S double a,b,eps;) K# @& Z7 D% ^9 I- D+ o
 { extern double frombf();
3 b9 Z: H6 p7 Y$ S; p8 g# p& K int m,n,i,k;
) I* e+ v$ N" i7 n0 @) q; G: t, [ double y[10],h,ep,p,x,s,q; b0 C: t+ h/ y5 Q0 L S. `
h=b-a;& s, ^9 U7 k2 M. k' r6 t
y[0]=h*(frombf(a)+frombf(b))/2.0;
% `7 ^) U! G% O# p# Y m=1; n=1; ep=eps+1.0;
) p6 l4 p+ }% G, |7 _8 ?6 n while ((ep>=eps)&&(m<=9))
' m9 K- N$ i/ j- W+ a { p=0.0;
+ n- w" y/ S; C( l for (i=0;i<=n-1;i++): A9 {9 T: V b, n0 S3 N% {
 { x=a+(i+0.5)*h;8 e7 t! H/ n8 R. K& P: m
 p=p+frombf(x);
) F# t4 r, ]8 A3 e# |( M }) L$ w+ |3 B; `3 }' O- y e$ N
 p=(y[0]+h*p)/2.0;' u- p/ H: C7 w& J% F" B4 E
 s=1.0;
9 t0 q9 k7 j0 j! p for (k=1;k<=m;k++)0 X( A% A: D9 |- l; d8 F8 o
 { s=4.0*s;
& g: ~8 {# f0 Z$ a( { ]7 X q=(s*p-y[k-1])/(s-1.0);3 {. Z* K5 `6 i/ d
 y[k-1]=p; p=q;# e0 X4 S) X9 Z. h3 O. n5 b# z
 }9 G) y9 Q7 L. P! F9 x
 ep=fabs(q-y[m-1]);5 H. b/ Z+ Q# L; E# ?3 F# ^
 m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
, p$ t" q3 z: U+ ?6 X }
, Q4 s1 a: [3 a4 `" ? return(q); D. ^8 r, N- s1 Q4 x: p) L9 {( E
 }<>呵呵希望对你有用！！