查看: 6001|回复: 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;9 ^$ |% [9 ]5 y! Q; ^
 }9 l9 v7 \+ Y8 h/ I, r% m
 }
' V0 k n6 P* K+ t0 [3 M. s u=(n-1)*n+n-1;. M4 U0 g( q! }* i& g9 H: N* x' ~
d=ar*ar+ai*ai;
- J8 ?9 H$ C, [1 \" [! _, r6 Y if (d+1.0==1.0)
6 z( w: u4 Y! L8 |) a: L2 x8 z) J { free(js); printf("err**fail\n");# s1 d2 C# I; W) y7 Y
 return(0);9 Z6 A, S# ~6 v8 d; a# I2 |$ `
 }$ M. }) B6 D9 L
p=ar*br[n-1]; q=-ai*bi[n-1];8 G+ y- ?- G$ }
s=(ar-ai)*(br[n-1]+bi[n-1]);
/ |, C2 l- G( }0 ]7 n br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
5 |: ?; x/ \/ g3 m8 i for (i=n-2;i>=0;i--)5 [! e# U- [. F# f! k) i
for (j=i+1;j<=n-1;j++)
! I2 Z; ~# T! @" R1 f6 x { u=i*n+j;
) N. `4 ^+ W+ ` p=ar*br[j]; q=ai*bi[j];
6 B0 m3 }8 h9 z o6 ] s=(ar+ai)*(br[j]+bi[j]);& i8 |7 Q9 D% N- B* O4 O" y: `# J
 br=br-p+q;& s3 t4 [! [8 c+ E( L7 n. D
 bi=bi-s+p+q;; z. s' t0 O6 A1 Z! L* t; b+ ~
 }
0 A, I9 h# `8 [! n( P js[n-1]=n-1;+ f6 f9 n5 t5 a2 P8 W% B
for (k=n-1;k>=0;k--)4 i- R s( b5 _- e! D( a3 z# b* D
 if (js[k]!=k)
2 p: H8 a6 V0 ] { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;- T' \! R" `9 t E0 H
 p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
 B5 A% H! V5 [ }
. L! ~, \; w, q' O$ w3 l' C4 c free(js);0 S: c2 G& ^$ B7 ~* _
return(1);- {3 ]0 g0 X7 x* A1 `4 G
 }<>平方根法<>#include "math.h"+ i% X+ h3 B6 ~$ l9 K- w% R/ I8 b
 #include "stdio.h"
9 ]. X- a( |# E/ a8 T int achol(a,n,m,d)
) p6 Y! x0 A0 h- }7 k1 B* w int n,m;0 i# J# n( K: ?
 double a[],d[];
, l+ P. M! S3 @+ x1 C& t { int i,j,k,u,v;
. Q3 x; J" `; T. a4 h if ((a[0]+1.0==1.0)||(a[0]<0.0))5 O- ~5 i Z1 v7 J. `3 G5 [
 { printf("fail\n"); return(-2);}
6 C S z G& ]0 b9 u2 K& }: S a[0]=sqrt(a[0]);# ^0 t+ C2 p. |
for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
, `& Y$ D1 L9 h for (i=1; i<=n-1; i++)
+ L" i% A p5 F$ P" _$ ^, G { u=i*n+i;
3 Z1 b! V% e3 z4 ^- o for (j=1; j<=i; j++)
; [& B; z2 E% a { v=(j-1)*n+i;, }! A c! s" C {, e% ~
 a=a-a[v]*a[v];
& Q. |' o3 G/ l& t3 Z }
7 [5 a# T2 U9 L( ~+ g; ~9 _( w if ((a+1.0==1.0)||(a<0.0))7 f) e9 c. k, e6 Z
 { printf("fail\n"); return(-2);}0 K$ ]: e3 Y# M5 h# F. L7 X
 a=sqrt(a);% s/ `2 x# E( [6 D; i+ S
 if (i!=(n-1))6 U( B+ l) `5 M0 ^" m/ M; l
 { for (j=i+1; j<=n-1; j++)
+ m# @: f* q1 M D# s. |2 B { v=i*n+j;
, y8 a# C* W" ? for (k=1; k<=i; k++)2 ?* S D2 P x; R
 a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];: r, W. ^7 J9 g
 a[v]=a[v]/a;
8 i, m" v v$ f/ J5 p }+ o7 w$ r8 P0 I$ n! A6 b
 }
; W" F+ z/ D1 g }
# y9 Y. w4 \- O$ ]' ? for (j=0; j<=m-1; j++)$ {1 c9 R- p6 E; L8 d
 { d[j]=d[j]/a[0];7 ?* b6 R% E* Y& O
 for (i=1; i<=n-1; i++); }7 R) x' c. G4 B9 }) @
 { u=i*n+i; v=i*m+j;! K* D3 V7 E) R2 M
 for (k=1; k<=i; k++)
) N4 ?! ]8 C+ r1 V7 p1 {$ q0 n9 I d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];; F/ c+ D7 `, @: a
 d[v]=d[v]/a;: T0 O8 a4 o3 b8 J4 i# n! X3 {
 }
' o" k* H/ n+ w* o' I% J6 A! V }9 B( [ [+ D' O
for (j=0; j<=m-1; j++)
& j+ d9 @! G1 ] { u=(n-1)*m+j;0 H. E* N1 B% ~- Y$ E% x0 L
 d=d/a[n*n-1];9 ? d* _3 g8 k& V
 for (k=n-1; k>=1; k--), h" H+ L- A" S; `; v+ N* Y
 { u=(k-1)*m+j;+ s, b4 B( s9 s/ T
 for (i=k; i<=n-1; i++): T" t" s! }1 m& Z/ d, s
 { v=(k-1)*n+i;; c) ]% K( v$ V% S8 g3 O6 Y2 i: \9 X
 d=d-a[v]*d[i*m+j];" k: S( U0 l- q7 M# l8 }
 }
3 @' J X+ Y2 K8 K# ?. y& p( S E v=(k-1)*n+k-1;
1 }2 b$ j* W$ i; N8 E/ v d=d/a[v];
1 Y$ d) B8 F5 d/ L' b) v } }$ |+ h6 z0 D2 g# G! f( N% n
 }
( P7 y4 m6 x% U" l* q' t" t: h return(2);
 K# p, s: h( u1 j# @% H" i }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t): y5 f* }7 @# l0 t1 W% Y
 int n;8 W9 u2 a$ G3 x6 U* B# i8 d. ?
 double x0,h,t,y[];
( K8 r. \ J4 n/ y2 ] { int i,j,k,m;
- U3 T; S, X+ B- T; u double z,s,xi,xj;
7 b8 g& b l7 c1 L+ ? float p,q;
4 I9 t8 _* w" }+ ` z=0.0;! A3 K/ ]4 x1 H% Y( m' A
if (n<1) return(z);! w; E _. N+ f2 @
if (n==1) { z=y[0]; return(z);}
+ r& S/ D# ]/ I5 T if (n==2)
; u, f& X2 P; k+ W { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;0 g" e$ H! x9 Z' p, Z; C9 Q
 return(z);3 M0 ?' C+ H" b% Z! `
 }1 h/ \' u, w9 R {
if (t>x0). k$ d6 Z% ?* H s$ |1 H ]" Q
 { p=(t-x0)/h; i=(int)p; q=(float)i;) e8 U$ b' L( d6 s5 R: Z
 if (p>q) i=i+1;! t4 v ^& X2 r! r( `# t$ ^
 }
: {3 ~9 q6 j9 w; j p) Z; g else i=0;
# X- `5 W6 V- p9 g k=i-4;3 k6 V1 D g' j$ e& |1 }, i8 d# U
if (k<0) k=0;
& m' F4 |9 h: Z m=i+3;: a/ y: X: K) o; Y4 B1 h
if (m>n-1) m=n-1;& {/ Y; X" M# a% G( E
for (i=k;i<=m;i++)! w% y# Q; A) c( m3 ?
 { s=1.0; xi=x0+i*h;5 s: i7 a1 Z* k' C; s5 V; F3 L
 for (j=k; j<=m; j++)' b! E6 k* r$ d. C7 P
 if (j!=i)0 i3 h P$ i0 ]2 Y% |4 B
 { xj=x0+j*h;
! X7 X" \" Q2 @& q& K4 Z, K4 ^( h7 H s=s*(t-xj)/(xi-xj);
) E$ ^1 N, K' I, G }% ?& h7 ]: A9 X3 V" j, T
 z=z+s*y;
0 }: t2 ]8 x" U# B8 g) ] }1 ]9 g' W& e# W) m5 ]7 u( P
return(z);
, A( N; `$ {+ L$ |( j0 | }) E7 \; ]' m4 W2 }, ?
向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"- d" ]2 M/ v: R/ w
 void hpir1(x,y,n,a,m,dt)
! b/ j/ D% c' |' l4 o8 L int n,m;
; `4 d4 X n* x1 T( ~3 r; k double x[],y[],a[],dt[];1 ?/ l, G* O2 K
 { int i,j,k;$ n& A: G! g! Y# B* m' m) y, `/ o/ |
double z,p,c,g,q,d1,d2,s[20],t[20],b[20];7 d9 r: T+ r, p6 E
for (i=0; i<=m-1; i++) a=0.0;
 A2 N2 b3 S8 z1 h1 ~8 j) U if (m>n) m=n;
, }. b8 o9 Z. O% V# O$ T' L if (m>20) m=20;
; C. j$ t& g; c5 i f8 J z=0.0;
* J! n, u4 k4 ]. P B for (i=0; i<=n-1; i++) z=z+x/(1.0*n);- k% ?5 C) I: z( O. c. I' f: A
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
8 k5 e2 @( ]" _7 u! z& o9 P V0 ^ for (i=0; i<=n-1; i++)
: H- A* o9 s9 w3 ]$ m' I7 O { p=p+(x-z); c=c+y;}# O; @) d ~- _! X/ O
c=c/d1; p=p/d1;. E7 }: K$ T( ^7 x
a[0]=c*b[0];$ A& m3 R) i5 W* `+ T% q5 W
if (m>1)
# V z! H+ C T' D { t[1]=1.0; t[0]=-p;' \9 D9 g: e, W( k
 d2=0.0; c=0.0; g=0.0;3 G0 E- v$ Y, Z i
 for (i=0; i<=n-1; i++)
& h% [ ?/ V4 V) u/ ?" U- ^ { q=x-z-p; d2=d2+q*q;- Z" j" a \* p) M% C4 e$ h& \% Y5 o% G; O
 c=c+y*q;
+ r( l$ D) W4 r- ?3 z' C g=g+(x-z)*q*q;- e" |% G0 B6 G+ d x
 }. Q9 Z5 ~4 M& B
 c=c/d2; p=g/d2; q=d2/d1;
! q, i/ O S/ r6 k! ] d1=d2;$ T) `" ^; r! J* z# V
 a[1]=c*t[1]; a[0]=c*t[0]+a[0];0 _. {) v: ^0 O; Z" ^6 s7 \/ X
 }
4 O$ Z, G. g8 r4 Y) f0 V for (j=2; j<=m-1; j++): C* q+ C0 [4 g6 @0 P7 B
 { s[j]=t[j-1];
) S3 c$ W h, l' H s[j-1]=-p*t[j-1]+t[j-2];
- {% r+ T0 [' ^+ s4 j0 L if (j>=3)
2 p; A0 g, J1 V2 B- z4 W( a5 b for (k=j-2; k>=1; k--)6 I$ k* ^$ Z0 \* S
 s[k]=-p*t[k]+t[k-1]-q*b[k];' A% z* U M' s$ s5 m0 c0 z
 s[0]=-p*t[0]-q*b[0];
 ?& ?: u3 ^% _+ ?1 q2 p$ D7 a d2=0.0; c=0.0; g=0.0;. |0 c3 R) S4 |
 for (i=0; i<=n-1; i++)# _- {# O: z& a/ g
 { q=s[j];4 F7 Z0 ^& J( q. B
 for (k=j-1; k>=0; k--)" d- D4 D5 }6 M0 f. \/ q! z
 q=q*(x-z)+s[k];/ J4 Z, w9 F I" t3 @+ a5 l W
 d2=d2+q*q; c=c+y*q;
( O2 l& n5 l5 ?! \ g=g+(x-z)*q*q;$ y5 m6 s3 A; _5 O5 J7 @4 K- V& o
 }8 ?6 e1 F8 ?& Q7 n% U. ~) G$ g. _2 r
 c=c/d2; p=g/d2; q=d2/d1;# h7 j: ?6 ]( @* l; L, E" w
 d1=d2;
9 u7 @- q( G8 G; v; ] a[j]=c*s[j]; t[j]=s[j];" s: r5 n1 O3 x$ T8 r' w. Q
 for (k=j-1; k>=0; k--)' v! ]2 k% t ~2 |
 { a[k]=c*s[k]+a[k];
! _# a0 h2 S5 M' B b[k]=t[k]; t[k]=s[k];
7 P0 g+ w* o$ h2 `) f$ Q }
5 F- M4 z( z w1 }6 C/ a' }. f }& a% F( \' G3 k. x7 F) L$ j9 \7 b- l
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
' s. K! J$ y' Z2 k* I0 N2 l for (i=0; i<=n-1; i++)5 `# j( _# P, i. r8 R& Z
 { q=a[m-1];
7 I8 f+ u. F/ y for (k=m-2; k>=0; k--)
3 m- H4 j; [% G; i q=a[k]+q*(x-z);
8 f4 W( v) P e {7 I p=q-y;* V0 Y. A- x, A3 x
 if (fabs(p)>dt[2]) dt[2]=fabs(p);: A, P1 ]* t; C$ D
 dt[0]=dt[0]+p*p;
6 m! D4 T+ T$ j# i, a) s8 s dt[1]=dt[1]+fabs(p);% W1 ^. B$ Y$ V5 W9 b
 }5 i" g3 A* N: ?" `" Z
return;8 v6 h+ c( b- U' {- g7 |' z
 }<>龙贝格积分法<>#include "math.h"7 o1 ?& |$ W+ F8 k
 double fromb(a,b,eps)* N* f: v% I& v: j* b' m2 s6 ^9 c, o
 double a,b,eps; s: Y4 T$ f' m! h9 v
 { extern double frombf();8 q* F4 T- [3 h8 Q! {" U
int m,n,i,k;
" H/ ^, O, S4 M- f, v5 w double y[10],h,ep,p,x,s,q;& t+ ]7 w5 c, d5 E; H( S; D" f# u
h=b-a;
/ R* J" C% [) S& O y[0]=h*(frombf(a)+frombf(b))/2.0;; u: Z7 V7 S0 w1 t; p- f- c
m=1; n=1; ep=eps+1.0;
, d6 m5 m- E( L* E( u/ q7 n B while ((ep>=eps)&&(m<=9))0 s l4 k3 J! X& R* C j2 @
 { p=0.0;# ~8 ^/ O6 ^0 X4 F Y$ k
 for (i=0;i<=n-1;i++), y$ X7 n1 _# m- g/ w! E
 { x=a+(i+0.5)*h;4 b" [: e6 N/ p& h2 C/ _
 p=p+frombf(x);
2 j* N1 a7 L) S3 F/ f }% y1 c( b5 [3 L9 ?- C! e: J( u3 v1 J, t
 p=(y[0]+h*p)/2.0;! ?1 C2 n( @9 N; m% ?
 s=1.0;8 z; f% }8 P! ^/ ^' B. a% L
 for (k=1;k<=m;k++)1 {! J+ v0 W+ V- |
 { s=4.0*s;
2 o/ x4 {3 d; I- v, q5 }% M* ` q=(s*p-y[k-1])/(s-1.0);
4 n, h" C- v' s y[k-1]=p; p=q;
! v+ p \" K$ i/ G5 t9 M }
7 R4 k' k# i0 r- m! X- h( f4 I8 B6 \ ep=fabs(q-y[m-1]);
; j2 m- O. y6 i6 }7 f+ f* w5 } m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
( w& s# k& E1 w6 j }/ x. C2 t( _1 w9 O
return(q);
. u0 c1 E* A% I0 s }<>呵呵希望对你有用！！