查看: 6004|回复: 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;8 Z' ?7 s* O" C |7 w1 G1 [
 }
2 V/ r/ O7 x6 L( _$ k! {) I7 I# | }
, O" {+ R4 Q+ ~. j9 T* r( N$ m) k% Y u=(n-1)*n+n-1;
$ l9 E( g! G$ c( o# I$ m$ G9 g/ C6 F d=ar*ar+ai*ai;
, \ v) D6 s6 l( ]4 s if (d+1.0==1.0)- `: y# o. f9 p
 { free(js); printf("err**fail\n");
( z$ g1 M) N3 J' ^' D return(0);& i* x) l. j4 C/ a; J: m* p! k
 }5 l) j6 L4 s' ~2 |( x, [" O
p=ar*br[n-1]; q=-ai*bi[n-1];8 ~3 L- T6 _' u9 m5 f& [
s=(ar-ai)*(br[n-1]+bi[n-1]);9 g8 L8 M5 b( B* R$ L$ a7 L9 H
br[n-1]=(p-q)/d; bi[n-1]=(s-p-q)/d;
; Q3 n' E# I$ s6 [4 l for (i=n-2;i>=0;i--) {) T& g4 n) A' [
for (j=i+1;j<=n-1;j++)
4 l1 }: O+ o Q1 ? { u=i*n+j;
 M1 F1 ~7 R5 J p=ar*br[j]; q=ai*bi[j];- {6 F X3 Q5 \9 Q5 _
 s=(ar+ai)*(br[j]+bi[j]);5 u# v: k$ _1 S( ~
 br=br-p+q;: r3 k) X( _2 R
 bi=bi-s+p+q;
* {9 s3 m) T6 d* O' d N$ m4 d }: F2 t- R' R* {% d. ^4 ^: ]
js[n-1]=n-1;
4 B: N; y" z* I+ v' D @ for (k=n-1;k>=0;k--)
; T2 @2 i5 n$ }. n3 l6 J! ~ if (js[k]!=k)
; L8 l* I: a( N$ D { p=br[k]; br[k]=br[js[k]]; br[js[k]]=p;
2 x; o3 ]* o0 _' b0 W0 R5 c p=bi[k]; bi[k]=bi[js[k]]; bi[js[k]]=p;# q" m7 q1 o0 ?: s
 }2 t; j# W% i) f8 a1 d) u+ w
free(js);5 l% B. d( o, r# t3 G
return(1);
 }<>平方根法<>#include "math.h"
) h) m) x5 T4 i) H) t #include "stdio.h"1 L$ W9 C+ I0 U0 \# U F. M- o
 int achol(a,n,m,d)
9 f. b# E+ t( S: g8 w2 n0 ` int n,m;
' e$ U& B8 P0 c2 a' A" N* F3 s" T! z double a[],d[];
/ h- v* {+ Y3 t! Z0 V { int i,j,k,u,v;0 ]$ N- v- @# h( W
if ((a[0]+1.0==1.0)||(a[0]<0.0))9 k& m' c2 @! j0 U4 g* K8 ^
 { printf("fail\n"); return(-2);}+ a8 A/ b! f, d M3 A
a[0]=sqrt(a[0]);
/ g, O% ~8 l1 y for (j=1; j<=n-1; j++) a[j]=a[j]/a[0];: Y3 @( R7 Q3 f% `" z
for (i=1; i<=n-1; i++)
# a3 U1 ?% G4 `7 z7 W- o5 d9 y" C { u=i*n+i;
2 i( l' x* q6 d2 i9 L# L& W$ [% A for (j=1; j<=i; j++)6 j6 w# [+ N- k1 f. ?2 v
 { v=(j-1)*n+i;% Y4 t8 ^7 B- C9 c. s4 c
 a=a-a[v]*a[v];: D* G3 I" w% P9 x5 J9 O1 @9 v$ I
 }
( A4 k- K+ a9 [4 K- F7 v% f if ((a+1.0==1.0)||(a<0.0)); j1 g! h( k; e9 N: ^" h7 \* B
 { printf("fail\n"); return(-2);}- g: [1 F3 u9 m0 v+ u
 a=sqrt(a);( V' p2 t* _! s- t2 e M
 if (i!=(n-1))6 U2 J& p/ P, O/ P; G
 { for (j=i+1; j<=n-1; j++): O1 h* a. }1 y; V6 e" Y1 n% D
 { v=i*n+j;
, H% w6 E9 P7 o0 t" i for (k=1; k<=i; k++)
6 f9 d) q: A* b. N4 X6 `6 X/ {0 ]2 O a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];& j$ D N0 x I: y# h+ |* P& o
 a[v]=a[v]/a;* I' p! l4 H' r% T* J0 Z
 }1 `7 W1 j2 i. P
 }
, A0 j1 g% d( n/ D0 y/ F }
1 v# g; a5 r( c; b0 t for (j=0; j<=m-1; j++)
# F& O; ^4 Y' P0 x5 M7 v { d[j]=d[j]/a[0];7 o7 t7 S# R* F3 n; K
 for (i=1; i<=n-1; i++)
) ` d. `, g$ A. v1 [8 |4 i { u=i*n+i; v=i*m+j;
! V9 M5 X, T% m" Q7 n _ for (k=1; k<=i; k++)5 D; V4 ?" q$ c
 d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
4 p" A' P1 ^2 C9 T d[v]=d[v]/a;1 N% R2 q# W1 x4 w% N
 }0 j' y" y3 b% k* C
 }0 j( N6 J9 O( `9 y$ N+ `4 ~7 q- l
for (j=0; j<=m-1; j++)- w, s9 m3 V! I
 { u=(n-1)*m+j;5 O* U5 n, o( H& g
 d=d/a[n*n-1];! s7 F' f- s9 Q7 {% a
 for (k=n-1; k>=1; k--)
. F1 w6 p3 m9 M7 [1 M, P& r { u=(k-1)*m+j;3 P3 M" B3 {/ R4 o5 ^+ s! C4 n! t/ ~
 for (i=k; i<=n-1; i++)4 c' E4 n" i5 V2 S5 J
 { v=(k-1)*n+i;
5 `$ L9 J; k T' _ d=d-a[v]*d[i*m+j];: S* `9 w5 ^" N$ Y/ h& i
 }) G7 X8 {7 m0 C
 v=(k-1)*n+k-1;/ q; t5 C; e) U( p0 a( y
 d=d/a[v];
: M& e# c0 U3 l; v/ N* ^ }
4 x7 h8 r% G7 I" T }" l. W3 @% V+ A& [$ L
return(2);
+ U$ v$ x; |1 Z9 p/ D+ B }<>牛顿向前，向后插值法<>double eelgr(x0,h,n,y,t)
3 V+ ?4 x: |( l, \# w int n;0 E7 @& D5 [& T8 w8 f, Y( e6 `7 R
 double x0,h,t,y[];
- I2 }8 j& ~$ Q8 l { int i,j,k,m;
( u1 n5 P/ _2 {' Q! s+ V double z,s,xi,xj;, Y' g! `' V; Q- D9 l" [
float p,q;/ `: A" ]$ H; Q+ J$ N! v+ \0 }
z=0.0;
. H$ y3 | v5 q* P9 @" [( x if (n<1) return(z);( n' L( j) l7 w# T5 Y: c4 M) u
if (n==1) { z=y[0]; return(z);}
; H2 i$ L) V( f2 L. k& j if (n==2)
0 ^+ ]1 _# ^( w" O) O { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;8 Z- ?7 v* O+ a) T. a; B1 S
 return(z);7 ], r' F% @% f+ D" v4 u" y) p* q
 }3 C" y7 y& M" v. b0 b0 L
if (t>x0)
+ e9 I' f1 x- I- K1 Q( w z* S { p=(t-x0)/h; i=(int)p; q=(float)i;5 j1 \9 Y% |3 g# u& ~+ w3 K
 if (p>q) i=i+1;% r) o+ Y) Q6 l9 H
 }1 s. a+ Q. ]# J7 P# p
else i=0;
+ \" m+ Q0 p3 D$ Z$ a2 h# s k=i-4;
+ w+ l% r7 P3 Q3 C9 |# O if (k<0) k=0;9 @% S0 b5 g% A
m=i+3;5 K. q: `) \3 w8 J
if (m>n-1) m=n-1;) u; ~; n. `5 O& [1 w, o+ U
for (i=k;i<=m;i++)
: a1 H& e! S! O* d { s=1.0; xi=x0+i*h;
/ o2 U/ T/ u; _8 {$ d for (j=k; j<=m; j++)
7 L6 x" l: L& i% ~8 i3 c0 J, c% @ if (j!=i)
' P/ t2 K1 U' p, H1 \ { xj=x0+j*h;
6 e& C$ @$ t8 L: c/ Z s=s*(t-xj)/(xi-xj);
! J6 A1 h& `: r! c! C/ f }
4 {2 d f6 x% { z=z+s*y;
/ ^9 m: I% w( o; \ z }) Z# w7 ~3 O, \/ V8 b1 F3 J
return(z);
# e2 t z. r4 ] `$ Y8 ~ }
6 G. [$ X: [6 P$ Q Z6 K向前，向后是一样的思想！！<>加权最小二乘法<>#include "math.h") j5 G0 r# g3 i1 K
 void hpir1(x,y,n,a,m,dt)
( h, I9 {& V7 ^4 F; l, o, s* I int n,m;
* E1 g9 Z! c2 W9 }# j% d$ r double x[],y[],a[],dt[];2 c8 ^1 [3 Y' e9 S3 B8 Q9 w& p3 I- E4 P
 { int i,j,k;/ Q) N/ \$ G3 l2 C
double z,p,c,g,q,d1,d2,s[20],t[20],b[20];
& c# N6 h) o) `3 u! W: i9 | for (i=0; i<=m-1; i++) a=0.0;
9 L7 ?4 @# P' o' E; X if (m>n) m=n; l" ?9 T* D/ t2 p r# S% e- ^
if (m>20) m=20;
, x, _. v1 g! `8 I z=0.0;- F! p+ o5 A$ {4 X! P6 m
for (i=0; i<=n-1; i++) z=z+x/(1.0*n);
/ `6 l' e* [9 a, k4 R5 P b[0]=1.0; d1=1.0*n; p=0.0; c=0.0;
" r" s7 |% O* F& ]9 G( T& l! [ for (i=0; i<=n-1; i++)1 [* n8 U' R- j% K4 g! ~
 { p=p+(x-z); c=c+y;}
5 k5 V9 G& i, n0 ~, e5 u# b$ s c=c/d1; p=p/d1;
0 G$ `0 k) `/ j2 ]- |' z* d. ? a[0]=c*b[0];. S& _$ T$ V0 g' p; Q& X) w
if (m>1)& w4 u+ V' l! @. G+ [- Q
 { t[1]=1.0; t[0]=-p;) E% L: p$ @, Q6 [$ P: _
 d2=0.0; c=0.0; g=0.0;
) I9 U& E, r- E/ F for (i=0; i<=n-1; i++)8 ] T+ X$ a6 E4 H* O
 { q=x-z-p; d2=d2+q*q;
0 F$ G# q( V C s1 p c=c+y*q;; D3 `* m4 E9 L2 ?- T
 g=g+(x-z)*q*q;
. K _6 }3 `2 \3 ]' g. r }2 {& V& F8 U# L4 J4 l' O
 c=c/d2; p=g/d2; q=d2/d1;% w" Z7 l m9 }+ U, V. p
 d1=d2;; c2 I Q; @3 z' q
 a[1]=c*t[1]; a[0]=c*t[0]+a[0];
/ N2 |( [& l- o$ W }
 z4 {% [3 ~! W5 K3 f/ n" z for (j=2; j<=m-1; j++). f1 p0 M' A0 j4 W5 [3 x! k
 { s[j]=t[j-1];
9 X- M4 |, r+ w3 T+ k$ L+ [ ~& \ s[j-1]=-p*t[j-1]+t[j-2];
! W& ]! T0 c5 e% q$ R2 m D if (j>=3)1 @+ @- ]1 P7 ?" T4 |: n3 ?9 {
 for (k=j-2; k>=1; k--)
. b& Y m( y: e- K. i s[k]=-p*t[k]+t[k-1]-q*b[k];! U7 r+ Y# `3 B
 s[0]=-p*t[0]-q*b[0];- d# X5 F4 D! M) \3 l* |
 d2=0.0; c=0.0; g=0.0;/ S: t4 C. R; g
 for (i=0; i<=n-1; i++)- _! p8 z1 S5 V, p' Z. h
 { q=s[j];
/ m2 m% I' L3 B! ?# G4 V: k: [ for (k=j-1; k>=0; k--)
& A9 Y1 D V( D' ~% q- a q=q*(x-z)+s[k];
2 ]6 i1 E7 I. m( K/ O& M d2=d2+q*q; c=c+y*q; A7 Y4 W6 |% k7 L
 g=g+(x-z)*q*q;
6 X) Z. U8 E* X4 |( J k% [+ g2 J }
$ }( t! F' @. C9 x6 c c=c/d2; p=g/d2; q=d2/d1; J3 X0 q5 H7 _: K3 }$ f i
 d1=d2;$ {2 ~* i/ J3 [0 F' r$ A
 a[j]=c*s[j]; t[j]=s[j];3 g0 F5 v' P" \( |. r7 c
 for (k=j-1; k>=0; k--)$ X2 L$ L _- ?4 e! o7 s# f
 { a[k]=c*s[k]+a[k];$ `1 V. m- C6 h0 p0 w; I* W
 b[k]=t[k]; t[k]=s[k];
4 E2 J4 ~( K9 b0 `; K( u }1 n6 u0 C- i8 }; b& q' F
 }
( |2 D. L8 @- ~9 S$ y dt[0]=0.0; dt[1]=0.0; dt[2]=0.0;
3 I. Q+ I2 g) L" S for (i=0; i<=n-1; i++)
1 t1 x" o6 @1 K4 ~9 [3 y4 x { q=a[m-1];
" N+ ^! t$ m5 h0 o- N7 B& \! @ for (k=m-2; k>=0; k--)
4 Z; i/ S& _! O8 P q=a[k]+q*(x-z);
' S2 V8 F) b. |9 r4 Q p=q-y;- \8 l- m8 U1 D
 if (fabs(p)>dt[2]) dt[2]=fabs(p);
) R: m6 J4 _. m/ w" ^1 L w* I dt[0]=dt[0]+p*p;) W# T/ F* v& S+ `; a* _
 dt[1]=dt[1]+fabs(p);
, G$ F5 W/ Z; B: E2 b9 x4 n: \& } }9 f- a# ?. u M, N! G
return;5 x2 r x( {- v+ u1 X9 {2 q
 }<>龙贝格积分法<>#include "math.h"
0 d& Y) d) \4 X4 c& M$ ~ double fromb(a,b,eps) Y% v2 K! E) e0 k
 double a,b,eps;
2 Y; A* e! b$ H8 X; d0 H @ `% C { extern double frombf();
* Z0 k* b% O9 h* S: {) T; ] int m,n,i,k;: [( f& Z5 R1 t7 {, T
double y[10],h,ep,p,x,s,q;, R& B" z/ J- b: D/ C4 a
h=b-a;
8 ^ m7 d* p9 @8 _ y[0]=h*(frombf(a)+frombf(b))/2.0;1 I1 b7 W3 ^$ `8 {; k& w2 k
m=1; n=1; ep=eps+1.0;( L4 W# e/ m- B# z4 m; o9 v
while ((ep>=eps)&&(m<=9)), ]1 K5 f O- f" A, x% d7 v
 { p=0.0;/ s$ u# X! c: L5 ~
 for (i=0;i<=n-1;i++)# `, O' `- ?, i( Z: x @" d
 { x=a+(i+0.5)*h; a* u+ N. s% u$ h
 p=p+frombf(x);
+ B0 |+ ]; I+ ]/ T# B% | }
- `1 h. B/ i9 K. O- ^ p=(y[0]+h*p)/2.0;6 q+ J+ a2 l* G& A2 D
 s=1.0;; `. M! O( L$ K! l1 O: d4 M
 for (k=1;k<=m;k++). f R1 E( P- o9 h/ U8 i
 { s=4.0*s;% S2 l, _( ~3 r/ s* a2 E$ O
 q=(s*p-y[k-1])/(s-1.0);% \$ a- f5 x0 b1 O) q$ L. V$ v* h
 y[k-1]=p; p=q;
+ X; `6 r2 R8 z. j" M: i8 L! k }( E( r+ K7 |" s* d2 w3 [
 ep=fabs(q-y[m-1]);
. I" q$ Q9 k6 v& | y) y8 k, W# X m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
2 C. T$ e, q5 ^& b+ x7 ?; Z }1 y* V; d* d( R, U8 H0 Y" L
return(q);
) O F6 S* O& S- {7 W }<>呵呵希望对你有用！！