查看: 5953|回复: 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;
) m0 V2 w7 c+ p( f }/ d4 I. d. f+ z# ]0 ?
 } J2 D( u7 z! ]( |) R: q
u=(n-1)*n+n-1;
7 N; x' `# z T0 M d=ar*ar+ai*ai;9 p6 c) {9 C8 T) ~7 C6 V
if (d+1.0==1.0)
4 u, X h( u* I- i0 h { free(js); printf("err**fail\n");
' z" y4 ~; Q; y$ P return(0);
0 [/ d' S8 A8 p: h% H& @9 T }
" s# L! M* ~( q" L7 A! q0 D p=ar*br[n-1]; q=-ai*bi[n-1];( p( h8 _/ `, k+ d9 x* G
s=(ar-ai)*(br[n-1]+bi[n-1]);9 Y; g% s( ^' d' X+ T/ G) K
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
/ P. I6 j: H! M for (i=n-2;i>=0;i--)+ z0 U6 E5 m4 E& |# e; c
for (j=i+1;j<=n-1;j++)
, L' S( U3 z; O3 @ { u=i*n+j;
8 T. r- b7 H9 l8 l( C8 A( |( T5 q p=ar*br[j]; q=ai*bi[j];
1 X; m: I N: B4 C s=(ar+ai)*(br[j]+bi[j]);
# E0 C- M& _* c* L0 g8 v br=br-p+q;3 T5 n9 e: _; h- C
 bi=bi-s+p+q;3 K& K1 z3 S2 C3 r0 `- |
 }9 G) W2 [8 B- O0 ]. q$ ^# {
js[n-1]=n-1;8 ` V5 R) A- Z4 ~7 y+ g+ [; ^
for (k=n-1;k>=0;k--)4 ?/ A8 T# |6 X3 F/ @
 if (js[k]!=k)2 {% U4 U' H6 Z) [
 { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;4 V/ P; b: e8 u5 x6 T; V K# h
 p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;) ?6 L) X) |0 F9 Q% g
 }3 f. a/ d2 X) D! K$ b: s; ?7 f3 A
free(js);
- H- U7 F9 c( {- Q return(1);
& Y, ~6 y4 Q8 r6 Y3 p }<>平方根法<>#include "math.h"9 B# M r. K- ]
 #include "stdio.h"
/ @4 w9 A# H+ h( B$ Y! J int achol(a,n,m,d)( S1 W5 Y, u$ k: p. T e
 int n,m;- I' T0 P8 \' u
 double a[],d[];, v: n9 \: t6 ]/ [6 @, B
 { int i,j,k,u,v;
$ u# ^ E' R* p) z+ t if ((a[0]+1.0==1.0)||(a[0]<0.0))
/ }; S+ ?. d' k5 ], A& L1 M5 z, c { printf("fail\n"); return(-2);}
9 g F) z6 G1 }. l* e a[0]=sqrt(a[0]);
; O6 [7 }" @4 H6 }+ v for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
: G% E' x4 l r. ^) \ for (i=1; i<=n-1; i++)
7 H. M/ f. {/ W: i { u=i*n+i;
) g7 q3 \! k$ b/ s for (j=1; j<=i; j++)
+ ]' q3 { H7 X) @: i { v=(j-1)*n+i;
 \( w) O/ X2 I" V2 I8 Q% ]) q a=a-a[v]*a[v];
0 @& |* q t a+ w }3 R/ k# F2 |) t; o' u- J& }: I
 if ((a+1.0==1.0)||(a<0.0)) v0 ~. S' Y% Q C. }
 { printf("fail\n"); return(-2);}
7 q5 R1 W7 ~) ?2 E+ W+ C a=sqrt(a);" V( n3 h; K: w+ P: s, G/ P7 p3 d% l
 if (i!=(n-1))
: l4 d1 \+ U6 c. ~0 k: Q { for (j=i+1; j<=n-1; j++)
( u; P W, S/ X0 j- e* M { v=i*n+j; k4 o. V; z4 f+ r! h% `, _8 h
 for (k=1; k<=i; k++)% w# A* m' Z3 _% A5 P) W" L3 N
 a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];" I4 g" D) K# z( v h H3 y0 Z
 a[v]=a[v]/a;' O% T, G/ `8 g: o& L
 }
$ C2 P( N+ g6 v' @' n' D }3 P+ i9 E, T- Q# W) R
 }
' T- F4 x4 D; }' [5 X9 o4 k for (j=0; j<=m-1; j++)4 V& f- }) @7 a( O' x6 J B
 { d[j]=d[j]/a[0];9 t: |, T `# v$ S% V: N% T) L
 for (i=1; i<=n-1; i++). ~" u7 C" }0 _( o" ?
 { u=i*n+i; v=i*m+j;) g& E% ^7 `1 E8 K
 for (k=1; k<=i; k++)
- _7 m/ K* O$ x# d5 l% m. K d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];% `, D b. J* \9 K) R* D E; {: j2 z
 d[v]=d[v]/a;
; ^& S8 q0 p6 `% F: L; x3 U }; c% N l% }. I6 N3 @6 ], Z: k8 O8 m
 }
, W5 k9 H- `; | for (j=0; j<=m-1; j++)4 n- G* M) I% G
 { u=(n-1)*m+j; ]0 }2 C7 f& h- {1 c% c E: E2 R
 d=d/a[n*n-1];) y8 T. v# Y' y( ^
 for (k=n-1; k>=1; k--)
( S/ {' a# X. H/ \/ J# B0 I { u=(k-1)*m+j;& j' P* ^3 _9 h1 h9 l8 |1 c$ k
 for (i=k; i<=n-1; i++)8 T$ R8 |8 r \0 Z$ Y/ T4 p
 { v=(k-1)*n+i;
7 K6 s6 `: b8 k, V d=d-a[v]*d[i*m+j];
" ]3 M1 m( a( B) {: r# ~ }" u7 U/ H) v1 @3 J. n0 {9 @
 v=(k-1)*n+k-1;
% j4 r# \" x4 ]$ h! n- Z d=d/a[v];
& b/ f3 L9 n/ j8 d }
$ d6 J$ G, Y4 e( Z( A }* t# z) Z: y- K
return(2);- l0 J' E& L! w+ A# E) z) G
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)5 Q; {+ }/ W4 V( w$ Y* R9 t
 int n;
# m8 p8 P; V R double x0,h,t,y[];% k* H7 l5 j6 D# |5 `' j
 { int i,j,k,m;% w N9 J6 d4 I: ^/ l5 \
double z,s,xi,xj;% G; t- V% l) N; }- J3 m/ h5 k4 L
float p,q;% ?* U; l" p6 q1 n& ^4 h0 N+ N. n* P
z=0.0;
! v9 i$ r8 m. w* H6 f( y* z if (n<1) return(z);
. I* N% w) w5 x0 k8 o5 Z( i4 m3 C if (n==1) { z=y[0]; return(z);}
" ]# U& T1 ^3 G if (n==2)8 E2 ^ P2 O: L7 ]7 A/ l$ ?
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
' K: M j0 H/ @2 M* ~; _# o return(z);
. P2 K4 [! {+ B6 K: H8 C }
9 Y6 U$ X9 o! S) n- J9 v+ @ H if (t>x0)
2 g3 e9 M* I+ K" d { p=(t-x0)/h; i=(int)p; q=(float)i;3 m' q2 h6 H% `3 W4 C' c
 if (p>q) i=i+1;5 q k* m3 p6 V# E
 }3 ^) T$ i0 U* |6 }& |3 q( Q3 ]
else i=0;7 E! F- B5 i9 H% o
k=i-4;, o4 S3 S- O7 H2 L
if (k<0) k=0;
' e; J G1 r2 y m=i+3;$ @6 \; {& D& G) N; q( Z7 N
if (m>n-1) m=n-1;
$ h* G* @% P! r' e for (i=k;i<=m;i++)# b6 v; K1 R% L
 { s=1.0; xi=x0+i*h;
+ @, q7 [ t0 n7 u: ^ for (j=k; j<=m; j++)
8 l+ a$ t, ~* ]1 V. m$ p2 Q6 t if (j!=i)9 ?7 d8 {9 f4 J
 { xj=x0+j*h;% i0 \; p( H8 Q5 x/ v
 s=s*(t-xj)/(xi-xj);6 n. c( Z9 i# ]2 U
 }0 l# T6 n; `2 o6 l: _8 |
 z=z+s*y;
: R* u4 U( |( n8 G' a1 z" K }
( T8 b* e9 ?3 M2 O" m" X- v return(z);
) |6 T$ O! M! z0 q# f, `/ z. ^ }/ f# ~, p! Z) d' @' [
向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h". b M ]4 b7 a9 P) P
 void hpir1(x,y,n,a,m,dt)8 A1 b" D* n3 I0 F& X
 int n,m;* T7 `$ |, v- d4 d# B- H5 u/ u1 o
 double x[],y[],a[],dt[];
5 I$ s, n0 H& X8 Q1 b { int i,j,k;
2 _# _' p) g) G9 O2 ~8 A double z,p,c,g,q,d1,d2,s[20],t[20],b[20];" ^+ H. k, w: e( X+ R/ U1 c
for (i=0; i<=m-1; i++) a=0.0;
* h3 s$ n5 @; Y/ d& |/ x! B; N if (m>n) m=n;
7 w+ m6 _# q8 P$ J. m7 o5 W5 E if (m>20) m=20; f8 [4 k. m z, n2 J! J) d
z=0.0;0 F k) s7 q1 `; w- O
for (i=0; i<=n-1; i++) z=z+x/(1.0*n);
/ T" A: j# X7 `: L9 W b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
: E8 u# C0 j% z for (i=0; i<=n-1; i++)/ b% A, S- k4 L8 p
 { p=p+(x-z); c=c+y;}; V( W2 |9 \6 r. T' [. S
c=c/d1; p=p/d1;: d7 e! J* Q1 q2 W7 u: Q: Q. p
a[0]=c*b[0];
; g/ S- r- p7 f+ k$ \: c if (m>1)
8 ^! c0 C$ R4 @6 S { t[1]=1.0; t[0]=-p;
( Q/ x7 {. W' u* i# j# d d2=0.0; c=0.0; g=0.0;/ u2 S! |( c9 I& E d/ G, r
 for (i=0; i<=n-1; i++)
% T! @+ \& {$ B/ z' | { q=x-z-p; d2=d2+q*q;
$ m# D0 v' j- j1 p c=c+y*q;- Z: M- z. w$ A
 g=g+(x-z)*q*q;: q% }7 ~% g" R
 }6 V. c/ w+ X1 u
 c=c/d2; p=g/d2; q=d2/d1;4 _ j+ F% Z8 I' X7 O+ N$ A
 d1=d2;
2 @( {9 ?2 N3 u' [* `! k7 G a[1]=c*t[1]; a[0]=c*t[0]+a[0];. D7 M$ j" j3 f& C& l C
 }' p# j. ~' \' R6 @) ?8 r3 u. e
for (j=2; j<=m-1; j++)' s& _; G' h3 } ^1 j3 [$ j5 u/ A. t
 { s[j]=t[j-1];, w# t; K$ t8 `0 v0 X" a G. p
 s[j-1]=-p*t[j-1]+t[j-2];
2 x8 c5 M1 ]2 t7 a" s if (j>=3)# F3 H$ R* ^& }
 for (k=j-2; k>=1; k--)
' A* S0 h* z7 O( r s[k]=-p*t[k]+t[k-1]-q*b[k];* v& z1 O& d g2 N: }
 s[0]=-p*t[0]-q*b[0];) M! Z/ F" t# J1 u
 d2=0.0; c=0.0; g=0.0;% N( w2 F2 V5 D4 q: J
 for (i=0; i<=n-1; i++)
) F8 d, a' }0 r- ]- I { q=s[j];4 G8 o! u( H7 k) ?# s
 for (k=j-1; k>=0; k--)
0 I- ?( N p7 ?) M q=q*(x-z)+s[k];
0 l' u; h- l- d9 x# l& D5 ^ d2=d2+q*q; c=c+y*q;
/ Z7 Z. \7 h" n4 n5 p g=g+(x-z)*q*q;
; M# p; k0 X* b. Y% E' C( D& w }
' \% p* D! R+ D! A% @+ s- }6 n$ |( W c=c/d2; p=g/d2; q=d2/d1;. v* a) T+ n3 b, o9 O( [
 d1=d2;
; [. X1 {+ ~. w a[j]=c*s[j]; t[j]=s[j];* s) u. f( H1 I# W
 for (k=j-1; k>=0; k--)6 n9 E, F+ C" s2 S" k' M# T! ]4 y
 { a[k]=c*s[k]+a[k];
 u) Q4 y5 `6 [2 l% B b[k]=t[k]; t[k]=s[k];; i# Q8 h% Y, Z4 r- f& `% i! Q
 }
6 N7 w) y. T" V& D# a7 y- H: j }
% E) T4 i" L/ m2 c& D- ~0 O1 T dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;( j) ?0 j7 Q+ a- M
for (i=0; i<=n-1; i++)
* `. x3 y/ P$ x: ]! x { q=a[m-1];
; V6 a0 T0 \4 }' E( E for (k=m-2; k>=0; k--)7 T! m7 {3 o. Y- I# c
 q=a[k]+q*(x-z);
# s" i- k+ F& p8 w x p=q-y;
+ q+ B5 n! G- {4 }% b if (fabs(p)>dt[2]) dt[2]=fabs(p);* B$ U$ P1 l! [: Q @, A
 dt[0]=dt[0]+p*p;0 E' Q3 F. `0 F& o1 d
 dt[1]=dt[1]+fabs(p);9 `2 q9 b/ a" }! m% y' T! M( ]
 }1 C6 K R7 }/ ?# F' ?% ~
return;
- h: N' [/ w& o2 ], N }<>龙贝格积分法<>#include "math.h"
5 W% ?0 j1 `* g# i: ?! X) E double fromb(a,b,eps)" U; B1 W2 i; B9 Z, m
 double a,b,eps;; T2 N) t: f! O3 s* p3 j& C& u
 { extern double frombf();
: t4 ^7 n; }5 x6 C( g& F int m,n,i,k;
0 E+ b$ X% M! `% o' x0 E double y[10],h,ep,p,x,s,q;+ m, M' o! h& g+ V% d
h=b-a;. a Y. u9 D; G
y[0]=h*(frombf(a)+frombf(b))/2.0;
' Y8 v2 s* I/ N/ \! B9 e m=1; n=1; ep=eps+1.0;. s+ ]/ v" t2 G% t G+ `: L
while ((ep>=eps)&&(m<=9))
0 o9 o0 S6 s. S5 s& M) \ { p=0.0;5 ?; D# R4 z. _
 for (i=0;i<=n-1;i++)
7 D Z% E* P/ i$ b7 g1 t { x=a+(i+0.5)*h;
# ]* T% b2 p A p=p+frombf(x);8 }0 J7 L" j. ^8 O" w$ k
 }
. h0 t( m0 j. c" K* r p=(y[0]+h*p)/2.0;! }; C, n5 k# w9 b
 s=1.0;
3 |, I0 t( ~4 e for (k=1;k<=m;k++): D l, Z9 k8 _* o; S% G
 { s=4.0*s;
+ r3 g8 W p# d: M" b q=(s*p-y[k-1])/(s-1.0);; n+ X# e: J7 t+ a/ B
 y[k-1]=p; p=q;% V3 [* b0 p( H8 q2 v4 T% v A
 }& I! n/ {1 Y6 t/ m! t
 ep=fabs(q-y[m-1]);
: S; d: n' ~8 Z) f# O l m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
% {9 Y* @ w2 }. A2 O. h/ d }8 V# {$ v. o( J% ^7 O3 s
return(q);
9 t# S3 H0 y2 Q9 I v4 p }<>呵呵希望对你有用！！