查看: 5962|回复: 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 K5 e5 o) k/ i }
; P: S% P! l0 p. q p+ P& ] }3 Y9 M2 y9 R% P' p X* u6 I ~ y
u=(n-1)*n+n-1;7 _' ?7 Q5 h2 C8 C# p
d=ar*ar+ai*ai;
1 e' k- ^* t: h/ y: S' y, g/ [# G if (d+1.0==1.0)# L4 [ N, \/ L. K
 { free(js); printf("err**fail\n");
+ X+ r- M" N- {, [$ J* N4 E& L return(0);& t/ \7 y) g& u; z3 l
 }
; ~! T- Q/ _) ?4 Q- t, i6 a# Z: P p=ar*br[n-1]; q=-ai*bi[n-1];
1 h6 a2 Q5 g7 }$ {; J4 A" V6 V! o. } s=(ar-ai)*(br[n-1]+bi[n-1]);4 m0 e1 @) Q- x0 d9 D# Y
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
. V) V7 z5 M9 x0 \: b- F. e for (i=n-2;i>=0;i--)
$ H7 }* p, n; n" V1 B1 b3 N \ K for (j=i+1;j<=n-1;j++)
) q7 c8 T9 }) ?0 B( o" E) [ { u=i*n+j;" L: P" a* y$ _) m( o2 ` k" V
 p=ar*br[j]; q=ai*bi[j];
. ^0 w- H) b7 M4 c s=(ar+ai)*(br[j]+bi[j]);
) h2 ~5 L; A4 f2 P | br=br-p+q;
/ Q6 h% ~ x5 C6 o. c" u bi=bi-s+p+q;6 c3 }9 ]" D' v) P$ Y! {
 }' g q; c& _7 H* V! M( b
js[n-1]=n-1;5 M) W* @6 h8 [8 d8 a* ~, C* J
for (k=n-1;k>=0;k--)
9 V$ y3 e, @6 u0 O* Z if (js[k]!=k)
 H8 }* c: R- \# e/ o1 N { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
3 F) E% L2 ] q; H p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;& s9 q. N$ F1 I4 I& `# ^
 }
 {) g3 e1 {8 M8 K free(js);
; q/ |, p; Q# g) y% F6 E/ C; z& M return(1);
/ g/ H% U# e6 x! W* n/ z }<>平方根法<>#include "math.h"
/ [; z4 o6 [5 ~2 y #include "stdio.h"5 m7 ~4 L$ H3 e/ M% S" [' C' b
 int achol(a,n,m,d)' `; O. C i* T3 n4 f
 int n,m;8 J# C: T( p( u3 J
 double a[],d[];% X9 ^, L8 \" y( _' C: K& X
 { int i,j,k,u,v;9 F: g1 Z \) L L: a" E# d
if ((a[0]+1.0==1.0)||(a[0]<0.0)), E [. l9 [% n( y5 q
 { printf("fail\n"); return(-2);}
6 n6 }6 E; K; M) D a[0]=sqrt(a[0]);
5 {# L0 a2 j1 A7 y6 }5 T0 S for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];2 U, f' e; {. E
for (i=1; i<=n-1; i++)4 I7 s3 J( ~' f) o! [5 N0 K: G$ K; S
 { u=i*n+i;4 T! y0 l$ s8 E0 z/ B
 for (j=1; j<=i; j++)
2 U; D# O ^5 ]: _7 ? { v=(j-1)*n+i;
4 [! e' J& k3 r# R" r- r a=a-a[v]*a[v];: w6 {7 X/ z& ]3 Q
 }, _, Q1 r6 H4 v" z, \
 if ((a+1.0==1.0)||(a<0.0))
7 M5 s, x" a. ]! H9 X# e { printf("fail\n"); return(-2);}
/ g& t' b R* Q1 ]+ Y% f$ y a=sqrt(a);5 m6 u$ z" S1 F; q r
 if (i!=(n-1))
9 |; D( r) h" K% G# L { for (j=i+1; j<=n-1; j++)- c5 L) P# y& U k; f2 q* A5 Q
 { v=i*n+j;
+ ?/ S% p+ b0 T for (k=1; k<=i; k++)8 d+ O- p/ [7 R) c: A# U6 h* Z
 a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];. f' b, M3 U' X5 v2 z6 S X1 f
 a[v]=a[v]/a;+ _( u* l8 u- r5 Y% j! e
 }. `* f' V6 M- l+ U i
 }0 w" }9 {- S% z- _* D2 N$ }
 }. l' Y( u" ]- J$ C0 ~$ F
for (j=0; j<=m-1; j++)
; A' q3 G; h7 O1 \: z: _ { d[j]=d[j]/a[0];
, H0 g0 G8 }% x- q- H1 Q+ B J for (i=1; i<=n-1; i++)* a4 a5 u$ @1 ~- j% Y
 { u=i*n+i; v=i*m+j;
0 m! O8 P. L- {" M9 p/ s/ O for (k=1; k<=i; k++)
; E2 t5 q6 c* _2 S; { d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];; t# f( X- Z2 e' B
 d[v]=d[v]/a;
" i8 D3 I3 I# a5 ^' R } L) B6 R! r/ A7 @& [
 }( J: j: e, A5 m
for (j=0; j<=m-1; j++)4 s! S: Q/ r2 C6 H/ l2 r: Q
 { u=(n-1)*m+j;
$ t/ k/ Z$ t0 l d=d/a[n*n-1];- m6 I R# s) L3 |) _
 for (k=n-1; k>=1; k--)0 [5 l+ Z3 w! W. d# O9 l7 p4 V( U {
 { u=(k-1)*m+j;' z6 F+ ~) t3 S) ?, ]$ D
 for (i=k; i<=n-1; i++)5 N1 C8 `1 Q/ ]
 { v=(k-1)*n+i;
. G- b: P- ?4 s. ~8 k+ @" z d=d-a[v]*d[i*m+j];- L- g9 [6 O+ ]5 D# R
 }
" ?# Z# k' b7 B0 c& v5 |" E! a, | v=(k-1)*n+k-1;- c( X0 @7 _& D; @/ f( A
 d=d/a[v];) e0 q/ f2 I1 M% q# e! g
 }
2 t4 C# u0 Y/ J1 A# C }- N D& @! U6 _( F: p* a P8 }5 {
return(2);
7 j, G+ V$ N- s }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)
5 P, E5 c7 z) c! B4 M int n;0 Q6 G& n; f9 k' E# ]5 h/ |2 d
 double x0,h,t,y[];
5 @0 ^% C" t. l; J# U7 b: r { int i,j,k,m;
4 F. Z- H3 V: k- L' Y6 R double z,s,xi,xj;+ b/ b' K, R6 E7 a" \. o
float p,q;# ~* M3 J& o8 C" M! L8 l
z=0.0;/ D. R% N( E# k
if (n<1) return(z);
 ` G ?2 x7 ~" A5 r* v8 U8 J: s if (n==1) { z=y[0]; return(z);}2 `4 p3 G1 |' [9 n: l
if (n==2)8 _, ~! z) L( B# d, k# C
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
- e0 [7 k6 z2 ~& C return(z);
9 J6 k, f1 o; o9 M# S7 i: Z }( _: J: o6 t6 l$ W5 @! b* T% s
if (t>x0)
. p9 j) p7 c5 P { p=(t-x0)/h; i=(int)p; q=(float)i;; M; N+ Y! V) n; a8 G' a. j
 if (p>q) i=i+1;' q z2 }& T3 A6 l
 }0 E, S; r9 F1 o1 w4 v
else i=0;
; \% I1 t9 M/ y: T k=i-4;7 g9 B) W( U. ]" m+ h& u
if (k<0) k=0;7 q0 }8 P) }/ \
m=i+3;7 G' u/ I7 O* h, E5 K5 C
if (m>n-1) m=n-1;
8 X. M9 @6 H- y) t+ [) u8 M for (i=k;i<=m;i++)
( s( ]' K& Q, a; c { s=1.0; xi=x0+i*h;; h" y, P ~' t
 for (j=k; j<=m; j++)% A' a! z9 v8 _, ]9 |2 R3 |
 if (j!=i)$ {6 q \, t) o; G9 p2 `
 { xj=x0+j*h;
' N, K0 c5 z g6 k( |& d s=s*(t-xj)/(xi-xj);, _) N* ~( G5 s5 ?* \4 M8 a' `3 g/ x
 }
9 q- F, _. y7 J z=z+s*y;
& \' @' k& R- N2 s4 S6 d% k! f }
8 O h, X# b* j) E% h return(z);0 A4 T( T( D4 V2 [
 }
4 S- v0 Q" ]; E1 k* c向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"
( R4 O* N, D7 Y; j, t void hpir1(x,y,n,a,m,dt)
6 U( s& E, o/ U% J6 |* o int n,m;" n& ]! L, ~' _
 double x[],y[],a[],dt[];
% ] a* V* `0 o+ C9 U { int i,j,k;) O |5 G9 ~/ j0 T# X# h* W2 g
double z,p,c,g,q,d1,d2,s[20],t[20],b[20];& s# n! i% I" f) d$ q+ K
for (i=0; i<=m-1; i++) a=0.0;" V' ]) F3 G2 D1 N" W0 L( e
if (m>n) m=n;
( p ^, S% o4 y: a8 q if (m>20) m=20;
$ w! }2 g x$ B% h+ N z=0.0;
; Z6 g9 [5 _% _ for (i=0; i<=n-1; i++) z=z+x/(1.0*n); {& B# y& Y4 g% f% o. H/ n7 V# ^+ K- [
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
" k9 Z3 k8 r/ J6 F* } for (i=0; i<=n-1; i++)5 P- h! }% a6 m* F: K0 z' ^
 { p=p+(x-z); c=c+y;}
- ~" t# O8 F! V6 c% g6 Q1 N+ P c=c/d1; p=p/d1;3 g5 M) @; \+ u( i
a[0]=c*b[0];
$ b1 y- i$ U5 M3 x4 d if (m>1)4 p0 t3 `% b4 V! I
 { t[1]=1.0; t[0]=-p;; b7 R% p3 y% M& _
 d2=0.0; c=0.0; g=0.0;9 h5 z% G6 m3 A9 t% v+ {
 for (i=0; i<=n-1; i++)
( _8 @ x6 v9 _0 t+ m { q=x-z-p; d2=d2+q*q;
, b# l' X+ Q' ^3 G, | c=c+y*q;
6 v; C+ H, d7 h1 b# B9 ]# ] g=g+(x-z)*q*q;
0 [) P9 O* R. o2 N" V5 b9 t* a( W }
 T3 X) F0 l% P+ T( o7 H! m c=c/d2; p=g/d2; q=d2/d1;6 R! }0 }5 D. }9 @
 d1=d2;
0 Z# W, T) t; { a[1]=c*t[1]; a[0]=c*t[0]+a[0];
( B5 x* y+ R. X' W/ j }
" s, i- a& |5 X0 V- h4 m for (j=2; j<=m-1; j++)
9 g# d0 d; z5 }+ B7 E { s[j]=t[j-1];1 j- x( g9 _% ]+ T% f
 s[j-1]=-p*t[j-1]+t[j-2];2 Q7 O& g/ E* l( M& W8 H% B% d
 if (j>=3)
5 N1 W) Y9 R& Q& Q for (k=j-2; k>=1; k--)
* `# [! z8 a* R4 N' N s[k]=-p*t[k]+t[k-1]-q*b[k];
! ?& k% M; w8 O3 k, C q s[0]=-p*t[0]-q*b[0];# ^; k! a2 D$ p. v2 {
 d2=0.0; c=0.0; g=0.0;! f4 u, y0 ~; C
 for (i=0; i<=n-1; i++)
- [: ?% v5 I- Q& h { q=s[j];
2 u, {2 p8 O: T for (k=j-1; k>=0; k--)
' ?) A5 \9 @! A' x# | q=q*(x-z)+s[k];+ e7 F v9 b H+ t9 o8 F/ O5 _- b
 d2=d2+q*q; c=c+y*q;
$ H, m7 @6 g/ |+ l3 H5 K! w2 u g=g+(x-z)*q*q;2 D- e4 J6 ?: c) M
 }+ [0 P5 q1 J1 y. X" o( n, m" @
 c=c/d2; p=g/d2; q=d2/d1;
; j4 w6 H6 h8 }2 f d1=d2;
% _" [/ J- Y3 A5 U( m3 z/ D a[j]=c*s[j]; t[j]=s[j];' N* r- S2 x. x# i$ j
 for (k=j-1; k>=0; k--)
' a8 c& r1 z" @ { a[k]=c*s[k]+a[k];2 S3 Q' Q' Z$ q B r8 T
 b[k]=t[k]; t[k]=s[k];* a# j! O% T) v" V' g. v4 r
 }% q5 h6 S4 [- S, Q" Q+ `2 G' x k! e
 }6 `7 o& c% e l
dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;0 K; S& Q0 S: T9 ^9 O" c: `- G
for (i=0; i<=n-1; i++)/ f5 x0 ?7 k" B, [
 { q=a[m-1];# j9 f8 {3 k5 C) a$ w7 K
 for (k=m-2; k>=0; k--). o/ T& H& U3 m1 I/ u
 q=a[k]+q*(x-z);
9 v9 d4 j; T! Y4 C# W5 i) r p=q-y;# c: O3 X$ Q5 B: z# U: r
 if (fabs(p)>dt[2]) dt[2]=fabs(p);' `: \( @& b: ?: {3 v' p$ f7 Z
 dt[0]=dt[0]+p*p;
- n- e6 H0 i* _! U$ y9 ^( h% k dt[1]=dt[1]+fabs(p);
4 x+ l% K4 t& t# F- V& X }, n& C* J6 a8 w: g; L
return;
: j0 B* b3 D% l: J8 E# o' o( z3 j }<>龙贝格积分法<>#include "math.h"* `3 R! |2 D, o( E
 double fromb(a,b,eps)5 q c" u3 r+ g
 double a,b,eps;- L. ^3 m% x2 i
 { extern double frombf(); ^* U6 v N* o" n$ x$ J% Q6 W
int m,n,i,k;0 i# \4 J$ h# U4 t
double y[10],h,ep,p,x,s,q;
3 s. ?8 `$ i+ ]1 _3 J( L h=b-a;+ F! g" a/ e+ ?1 l
y[0]=h*(frombf(a)+frombf(b))/2.0;. ^2 c! Y& g5 t* s) Z% v
m=1; n=1; ep=eps+1.0;% Z& F/ H/ A7 h+ [ k; w! B
while ((ep>=eps)&&(m<=9))6 W! b; p7 @& g$ b" k, U
 { p=0.0;
; p" D3 v! w2 b! x" N$ M& k8 i for (i=0;i<=n-1;i++): ?3 {& e+ H1 N9 n$ {! b! D
 { x=a+(i+0.5)*h;# l- D9 `9 F# [/ n1 x
 p=p+frombf(x);! `) s$ q$ K+ i7 |
 }* ^ ^" d; A( ]) ~2 `" r$ {1 U
 p=(y[0]+h*p)/2.0;
4 ?4 N, Z* d; n" G: q s=1.0;
' ]/ v/ P! z3 e for (k=1;k<=m;k++)* w. N7 D: ~7 r; O6 P y/ K
 { s=4.0*s;
, K( t9 {7 _0 W* U- l# ~% N q=(s*p-y[k-1])/(s-1.0);$ t2 S: w. I6 ~: j; o" L3 \
 y[k-1]=p; p=q;
' B7 m: A& n! |+ c: Q }
' l; A4 E. N/ h& p. x ep=fabs(q-y[m-1]);7 `+ U P$ X7 l" ~5 }/ O
 m=m+1; y[m-1]=q; n=n+n; h=h/2.0;7 ~. f# d5 v3 }! S
 }; y- R9 Q" @. M r# j' ^: T0 Q* N N7 d
return(q);7 L7 u2 E, j' `- r: {
 }<>呵呵希望对你有用！！