查看: 6039|回复: 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; @" ~& V w0 I9 I$ ~1 i! [
 }
' C0 S3 O+ H( P8 c! s+ ?- t }
' u0 a3 [5 q C- X7 j u=(n-1)*n+n-1;/ I, f) x0 F3 z$ M; S- ~
d=ar*ar+ai*ai;2 |3 E J9 S/ u( ]
if (d+1.0==1.0)
) z* `" J7 n5 ^+ r { free(js); printf("err**fail\n");2 ~; ?" U. @$ O0 C2 m6 F$ o
 return(0);
4 |$ E: a6 F) v$ Q; \ }& E H& N9 ^- C+ ?: r
p=ar*br[n-1]; q=-ai*bi[n-1];
6 [/ R/ F2 b9 @8 S- O# R s=(ar-ai)*(br[n-1]+bi[n-1]);+ \0 A }# k* w* U- m
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
/ W2 [ }% h+ Y3 K6 m for (i=n-2;i>=0;i--)
6 ?. b H0 s( S. d% Z o* C for (j=i+1;j<=n-1;j++)
4 d& v" L, @; ~3 H" U3 ?. ~ { u=i*n+j;! s6 u1 X( f. Q# g
 p=ar*br[j]; q=ai*bi[j];
, n$ y( p; O- b9 N' ~ g/ S/ d2 p0 l s=(ar+ai)*(br[j]+bi[j]);
5 d) G6 z% n8 |6 `8 A% A U; C2 k br=br-p+q;! d3 V7 G5 Y, I- l
 bi=bi-s+p+q;: b" {( `1 A& _6 B) I
 }
$ v7 X+ G( `% {& S0 S js[n-1]=n-1;4 j- C6 I% A; c& x2 _. W
for (k=n-1;k>=0;k--)1 O6 b. H9 W# f) z6 |% b
 if (js[k]!=k)
6 f" l" a& A* w; v { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
. ^. o* Q6 N* t3 o/ y% V3 s p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;
3 x4 y4 m @! {2 X) U* H }
; {+ }* z% c$ Z- P( e7 w* E free(js);8 U+ `, u) G: |2 L: S& w! i+ N
return(1);' {3 [/ A; d7 @- O
 }<>平方根法<>#include "math.h". E8 z) U4 M! f8 C' X* P9 j j
 #include "stdio.h"
/ f1 }, x/ c: u/ u4 B" _ int achol(a,n,m,d)3 [- S$ _. A! \% W' B
 int n,m;
0 i- [8 J, V) ?, o/ E8 [, Y double a[],d[];
9 K( e, r8 M" o { int i,j,k,u,v;
$ I+ z/ C( V4 R$ v. T if ((a[0]+1.0==1.0)||(a[0]<0.0))' S1 b$ x* S1 q- ?0 w
 { printf("fail\n"); return(-2);}. p8 G8 r1 V" u$ T/ p
a[0]=sqrt(a[0]);; s8 F$ M F! k. k
for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];
3 T& k. Z; y" f% o for (i=1; i<=n-1; i++)7 h" U5 H- @1 P! w$ Y2 t
 { u=i*n+i;
; k" e1 s2 t+ x2 }# u( @ for (j=1; j<=i; j++)
" l7 @) h# Z3 K' k, A { v=(j-1)*n+i;# d1 k) j' ?* Q& g2 |) O
 a=a-a[v]*a[v];7 Q! t- l C* W1 F @
 }4 x& ~$ s" i1 U! x9 P
 if ((a+1.0==1.0)||(a<0.0))! X. P5 N. \3 G
 { printf("fail\n"); return(-2);}- c* M5 w+ X8 V
 a=sqrt(a);
7 S! m$ \8 Y) i9 y; C" d3 |, ] if (i!=(n-1)): H$ J7 U! V& \( h8 s+ J
 { for (j=i+1; j<=n-1; j++)$ L1 q5 c( }3 W4 s, L5 D
 { v=i*n+j;
! C7 g6 C$ d9 Y( R for (k=1; k<=i; k++)
3 U( e, B ^; E a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
5 P0 F, e, q! O( w) i) _8 O) t; R/ m a[v]=a[v]/a;! ]9 n4 ~- f' k5 \( }" w
 }
) ]! B5 R- J' ] }
7 c6 m; ]& E: {$ o8 B, ?1 O2 z }
# h ~: W% [+ V" ~9 ] for (j=0; j<=m-1; j++). t. H! u6 F) D2 c) K' r2 u( R, o
 { d[j]=d[j]/a[0];. V1 `. F$ S0 t$ F( K
 for (i=1; i<=n-1; i++)0 E C# @; `6 x) ^% c3 o& ^. o
 { u=i*n+i; v=i*m+j;) M8 c3 D2 ] B5 s/ J/ s
 for (k=1; k<=i; k++)' Z, l1 y6 l" }/ R/ @$ y% @: i6 Z' E! J
 d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
6 S/ @. d6 t8 o: E) ? d[v]=d[v]/a;
! g& O7 L) @" w' a9 `; j1 N }
# @+ }9 t) `# |6 M/ G9 H1 \ }6 R$ f& X0 y% Q' e) {7 d4 Q
for (j=0; j<=m-1; j++)
/ J* t( N( t1 n7 B' _ { u=(n-1)*m+j;9 M6 w" \: }) O9 @7 h5 V7 x# ~
 d=d/a[n*n-1];& @0 p+ w, K0 c0 x y
 for (k=n-1; k>=1; k--)! z# j9 r" Y5 W4 N$ r/ K
 { u=(k-1)*m+j;
9 X9 O6 w+ K$ e6 V8 K for (i=k; i<=n-1; i++)3 _ r9 v0 u6 D. N7 c5 K% N
 { v=(k-1)*n+i;+ _; u" [8 V) B3 K: [" P1 C
 d=d-a[v]*d[i*m+j];+ h( y3 }1 t9 z" D- {4 [7 e
 }0 n: Y7 q% w- n- C1 W
 v=(k-1)*n+k-1;
; p. f" N3 e: f$ ~. W8 F4 ~ d=d/a[v];8 M6 g/ n4 w1 v6 C" q8 `, ]. |
 }
: \( e* ?7 d9 [/ r9 P }* p7 j8 Z6 q! t
return(2);$ J) K: R7 Z1 H3 B; s4 l
 }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t) f* Z$ m4 P/ r
 int n;8 Z: h6 K" B `; b) k9 [9 \
 double x0,h,t,y[];5 K; G* e1 O7 b+ g+ y7 o9 R
 { int i,j,k,m;
9 j0 O" R5 P' r1 y5 J double z,s,xi,xj;# n$ h8 }! u9 `1 ~* K
float p,q;( `$ ? i4 y5 s0 W {; G
z=0.0;
) d: Z3 S! S7 G u$ N if (n<1) return(z);2 R0 A. I: [, p1 e+ V% `
if (n==1) { z=y[0]; return(z);}. ?" l' _, x3 Z i. |1 [
if (n==2)% B& |- }, |8 \* k/ V
 { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
. U1 h2 X) S+ [. Z return(z);
, @ Q1 V9 h: h$ w }9 n. J/ d5 R1 j" }4 V: r
if (t>x0)( S3 _# ~( d5 N. _
 { p=(t-x0)/h; i=(int)p; q=(float)i;
& X! Q1 {6 n/ c+ I if (p>q) i=i+1;$ {6 T( {/ ], j; a. [+ w
 }1 r8 b- [3 Z- f
else i=0;
3 _6 s2 k8 B2 ^- |! P k=i-4;
7 E. V& p0 V/ Z' x: l z if (k<0) k=0;8 ?" g9 }+ ?# o! f: p$ k+ [
m=i+3;
9 w0 w7 H' a3 k5 Q J if (m>n-1) m=n-1;1 o: R. b/ ]1 O$ S1 i
for (i=k;i<=m;i++)
. Y2 H# I. C- ^. x* { { s=1.0; xi=x0+i*h;. ~& z- B) I" [, @. |
 for (j=k; j<=m; j++)
" V; a/ U) @ A* F if (j!=i)
$ o. n( G+ B3 x y { xj=x0+j*h;0 G8 z! U1 s( m' D- h: Q' ?! W
 s=s*(t-xj)/(xi-xj);7 Z5 K% g5 R2 J: C6 Q2 Z: _
 }
5 Q* A" G3 i. S4 I' P+ k! R z=z+s*y;: b9 g* k. Q y1 r4 B. u) D9 n* L
 }7 j z+ T2 }+ U
return(z);
+ Y( E& ]. Q( Y Q# F L }
* u- Q) b% K* r \* T向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h"8 [1 C9 @9 m: D$ p* z, W
 void hpir1(x,y,n,a,m,dt)$ h$ e) ^. P2 L, E* L
 int n,m;
: d& t) s; V3 V: ^5 E5 |, J" ` double x[],y[],a[],dt[];- M' U7 t' C, I% t
 { int i,j,k;
2 a3 G% Z1 _, T+ g2 u1 e double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
 j, C( I( h' A5 G4 m+ ? for (i=0; i<=m-1; i++) a=0.0;
/ b' G( b. I6 Y2 h8 V6 _* ] if (m>n) m=n;# r$ X j' a0 v3 A
if (m>20) m=20;& a5 p& G/ p3 n9 e1 E& i( ^
z=0.0;5 F/ H y" R# z) ?0 J$ R
for (i=0; i<=n-1; i++) z=z+x/(1.0*n);6 d: K8 U5 A, {5 w6 e
b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;7 l2 H7 p) f% b% `/ R
for (i=0; i<=n-1; i++)( ~) z) k% Z: T+ D8 G% H
 { p=p+(x-z); c=c+y;}
* l3 Q" M! G- V1 \ c=c/d1; p=p/d1;
5 Q0 s7 C4 z2 w a[0]=c*b[0];6 m" M% X7 X) {5 D' u0 h% m9 ?
if (m>1)' Z6 \% Z# @/ q, L$ }) ^( z/ t. t
 { t[1]=1.0; t[0]=-p;
9 C% O$ {- g) W/ J/ i+ m d2=0.0; c=0.0; g=0.0;
/ y6 K! C8 k/ U( N# R for (i=0; i<=n-1; i++)* z1 l1 N2 E+ y
 { q=x-z-p; d2=d2+q*q;
4 z+ Z9 ?, e( Q+ A( ?0 l- U, _- C c=c+y*q;% z: a1 c" `2 B# |7 e
 g=g+(x-z)*q*q;
0 l4 A, ?. h' l U; B }8 f/ x0 W0 ~/ n, [, h
 c=c/d2; p=g/d2; q=d2/d1;6 f8 \7 \+ ~0 L4 f: J
 d1=d2;# H. g' s8 B1 [
 a[1]=c*t[1]; a[0]=c*t[0]+a[0];$ h4 I5 {1 g5 E/ z. W
 }
' k3 W, T- F: N/ b+ ?) A! @1 L for (j=2; j<=m-1; j++)) d1 t1 a- b) A$ T9 p; V+ X
 { s[j]=t[j-1];
3 e5 ^# q& `8 Z4 [8 K1 D6 t s[j-1]=-p*t[j-1]+t[j-2];
) H4 j8 T2 g4 M if (j>=3)7 R: D+ R% W7 r3 X8 V q @
 for (k=j-2; k>=1; k--)% ]% {) ?' `6 ~. _3 ]+ E
 s[k]=-p*t[k]+t[k-1]-q*b[k];
& ~2 l- Y2 U1 f: i: `% @ s[0]=-p*t[0]-q*b[0];' w% ?) ?* K" |& e
 d2=0.0; c=0.0; g=0.0;# Y/ P0 {8 W4 T z/ F5 ^
 for (i=0; i<=n-1; i++)
' n; \. @1 ?/ J2 {- w3 O+ N { q=s[j];
* c3 L. C1 `& U. U8 { for (k=j-1; k>=0; k--)
% d- N7 Q& y) N" [3 ^2 f q=q*(x-z)+s[k];
 F; A4 l: W* i& F- a d2=d2+q*q; c=c+y*q;( M* a/ D7 Q. u; q9 G
 g=g+(x-z)*q*q;
* a) c( ~& E- b5 u/ V! J: v }
, ]. H! T, q& w0 M* S c=c/d2; p=g/d2; q=d2/d1;
, s; A- d5 G4 X: | d1=d2;* ~2 ]2 \& v0 a. w% z+ ]/ Q
 a[j]=c*s[j]; t[j]=s[j];
- P; a( ~5 g( X5 A# Z. ] for (k=j-1; k>=0; k--)
; g2 m/ b, J) O6 W. H$ N { a[k]=c*s[k]+a[k];5 ^+ G8 {# f3 n$ p' c' V, I
 b[k]=t[k]; t[k]=s[k];
/ ?2 @1 Y0 u! X0 r# F6 x& y }& J; R9 R, l1 j0 T" B* ?$ W3 E
 }
, j7 M$ j1 k& R. P1 l& z, j dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
4 N5 a$ j4 v4 r" b* k. m5 ]' d/ { m for (i=0; i<=n-1; i++)7 z4 g3 I8 t# D6 _4 S
 { q=a[m-1];
5 ^3 X! u1 H, d; j2 u) f: W5 i for (k=m-2; k>=0; k--)2 ^9 T* L1 I2 R0 D
 q=a[k]+q*(x-z);
' I8 B& d( v( D1 o8 h2 D; i p=q-y;
( R# ^* @6 w4 L' S) P if (fabs(p)>dt[2]) dt[2]=fabs(p);
' _2 \0 M8 P2 j# ^4 q* C( C dt[0]=dt[0]+p*p;- S9 \8 S3 @$ A y3 }$ e' h
 dt[1]=dt[1]+fabs(p);' O- ~# M& y g
 }
5 T1 l' `) l) L+ A$ E return; v8 r! M9 v* u5 @
 }<>龙贝格积分法<>#include "math.h"8 J7 V8 M. C8 k; v& j
 double fromb(a,b,eps)
' K5 o+ r! V* T# x0 ~4 A double a,b,eps;
* t8 h. t/ _& C, R { extern double frombf();. ^, T8 |& }: q2 G# P M2 s2 s) V
int m,n,i,k;7 \7 W6 R1 o/ L) ]6 F( {1 U6 B# f
double y[10],h,ep,p,x,s,q;
) e* B1 K5 m8 J$ x h=b-a;0 ~* q# U5 U5 n$ { B- D
y[0]=h*(frombf(a)+frombf(b))/2.0;- E9 @/ z3 H! x' u0 J
m=1; n=1; ep=eps+1.0;
9 b! e# p, f2 k- H+ X: o [2 s while ((ep>=eps)&&(m<=9))
! ~% K% T. L0 N$ A( C" y7 P { p=0.0;3 a, ]/ l0 W& n0 l
 for (i=0;i<=n-1;i++)
$ a" ~1 |. A5 d { x=a+(i+0.5)*h;
+ q2 P( h; {* U8 _% @ p=p+frombf(x);
- a) W4 R2 f- ^* M% P }
6 L# ]: u& \8 o- p0 { p=(y[0]+h*p)/2.0;5 ~& f" _/ X' _& r" K1 l
 s=1.0;. v# l$ E0 a8 R$ m
 for (k=1;k<=m;k++)( k+ W6 U j6 Z- [1 ~% x
 { s=4.0*s;7 z2 `- q# O1 U, v( |- `; t
 q=(s*p-y[k-1])/(s-1.0);, V% r; f1 K4 e3 r- z. f: ^
 y[k-1]=p; p=q;
& N& u6 `: I% ]/ u1 {) \ }
9 W% w% {& O- S& j7 b0 ^0 `2 } ep=fabs(q-y[m-1]);
$ l5 V( N. w6 y1 L* u4 B m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
1 Y+ E {% F( B& c# J }1 K5 n7 ~7 @) M: Z/ b/ B8 v. T
return(q);! G) e/ \1 u9 l. j0 \7 y! k( s
 }<>呵呵希望对你有用！！