平面矩形Gauss數(shù)值積分

1、平面矩形Gauss數(shù)值積分/計(jì)算函數(shù) f(x)=x*x*y+x*x+x+1+y*y+y; 在a,b上的積分 , 采用每個(gè)方向上3個(gè)積分點(diǎn)#define a 1#define b 4 #include "stdio.h"#include "math.h"double gaussx9,gaussy9,H3,gausspoint3;void main() int k,i,j; double result,fx,x,y;H0=0.555556; /積分系數(shù) H1=0.888889; H2=0.555556; gausspoint0=0.7745967; /積分點(diǎn)

2、坐標(biāo) gausspoint1=0.0; gausspoint2=-0.7745967; for (i=0;i<3;i+) for(j=0;j<3;j+) gaussx3*i+j=gausspointj; gaussy3*i+j=gausspoint i; result=0.0; for (k=0;k<9;k+) y=gaussyk; x=gaussxk; fx=x*x*y+x*x+x+1+y*y+y; result+=Hk/3*Hk%3*fx; result=(b-a)/2*result; printf("%f", temp);解線性方程組之 LU分解法#

3、define SIZE 4 /積分點(diǎn)數(shù)#include "stdio.h"#include "math.h"void lu(a,b,n,l,u,x)float aSIZE,bSIZE,lSIZE,uSIZE,xSIZE;int n; int i,j,k; float sumb,suml,sumu,ySIZE; / for(i=0;i<n;i+)/ u0 i=a0 i;/ l i0=a i0/a00;/ for(i=0;i<n;i+) for(j=i;j<n;j+) sumu=suml=0.0; for(k=0;k<i;k+) su

4、mu+=l ik*ukj; suml+=ljk*uk i; printf("%fn",l ik); u ij=a ij-sumu; if(i=j) l ij=1.0; else lj i=(aj i-suml)/u i i; uj i=0.0; l ij=0.0; /前推過程 sumb=0.0; for(k=0;k<i;k+) sumb+=l ik*yk; y i=b i-sumb; /回代過程 for(i=n-1;i>=0;i-) sumu=0.0; for(j=n-1;j>i;j-) sumu+=xj*u ij; x i=(y i-sumu)/u i

5、i; /主函數(shù)int main()static float aSIZESIZE=1,2,3,4,1,4,9,16,1,8,27,64,1,16,81,256; static float bSIZE=2,10,44,190; float lSIZESIZE,uSIZESIZE,ySIZE,xSIZE; int i,j,n=SIZE;printf("input N(N<30)");/ scanf("%d",&n);/* for(i=0;i<n;i+) b i=0; for(j=0;j<n;j+) if(i=j) a ij=2.0;

6、else if (fabs(i-j)=1) a ij=1.0; else a ij=0.0; b0=1.0;if(n/2!=0) bn-1=-1;else bn-1=1;*/lu(a,b,n,l,u,x);printf("The left matrix is:n");for (i=0;i<n;i+) for (j=0;j<n;j+) printf("%6.3f ",l ij); printf("n");printf("The ultra matrix is:n");for (i=0;i<n;i+)

7、 for (j=0;j<n;j+) printf("%6.3f ",u ij); printf("n");printf("The solution is:n");for (i=0;i<n;i+) printf("x%d = %6.3fn",i+1,x i);解線性方程組之 LDL'分解法/LDL'分解void LDL(double ANPNP,double XNP,double fNP)int i,j,k;double sum1=0,sum2,sum3,sum4;sum2=sum1;su

8、m3=sum1;sum4=sum1;for(i=0;i<NP;i+) for(j=0;j<NP;j+) L ij=0;for(i=0;i<NP;i+) d i=0;for(k=0;k<NP;k+) sum1=0; for(i=0;i<k;i+) sum1+=Lk i*Lk i*d i; dk=Akk-sum1; for(j=k;j<NP;j+) sum2=0.0; for(i=0;i<=j-1;i+) sum2+=Lj i*d i*Lk i; Ljk=(Ajk-sum2)/dk; for(i=0;i<NP;i+) sum3=0.0; for(j=

9、0;j<=i-1;j+) sum3+=L ij*Zj; Z i=(f i-sum3)/L i; for(i=0;i<NP;i+) Y i=Z i /d i;for(i=NP-1;i>=0;i-) sum4=0.0; for(j=i+1;j<NP;j+) sum4+=Lj i*Xj; X i=(Y i-sum4)/L i; for (i=0;i<NP;i+) printf("n U%d %7.4f",i+1,X i);解線性方程組之 選主元Gauss消去法void chg(double arr_aNPNP, double arr_bNP, int

10、 m) int i,j,k; double a_chg, b_chg, max; max=fabs(arr_amm); k=m; for(i=m+1; i<=NP-1; i+) if (fabs(arr_a im)>max)max=fabs(arr_a im); k=i; if(k!=m) for (j=m; j<=NP-1; j+)a_chg=arr_amj; arr_amj=arr_akj; arr_akj=a_chg; b_chg=arr_bm; arr_bm=arr_bk; arr_bk=b_chg; void gausszy(double ANPNP, double

11、 fNP) int i,j,k; double a_sum=0; for(k=0; k<=NP-1; k+) chg(A,f,k); for (i=k+1; i<=NP-1; i+) A ik=A ik/Akk; for(j=k+1; j<=NP-1; j+) A ij=A ij-A ik*Akj; f i=f i-A ik*fk; XNP-1=fNP-1/ANP-1NP-1; for(k=NP-2; k>=0; k-) for(j=NP-1; j>=k+1; j-) a_sum=a_sum+Akj*Xj; Xk=(fk-a_sum)/Akk; a_sum=0;

12、另外一個(gè)/ANPNP 是系數(shù)矩陣，fNP是右端向量。/分為消去和回代兩個(gè)過程。結(jié)果存放在XNP中，是解向量。應(yīng)為PUBLIC型/這里驗(yàn)證了LDL和gausszy兩種方法的差異，LDL要求嚴(yán)格對(duì)解占優(yōu)，/而Gausszy不要求。#define NP 3#define N 3#include "stdio.h"#include "math.h"double ANPNP;double fNP,func;double XNP,YNP,ZNP,LNPNP,dNP;#include "gauss.cpp"#include "gaussz

13、y.c"#include "ldl.c"void main() int i,j,k; double sum,sum1; extern void gausszy(); extern void gauss(); extern void chg(); extern void ldl();/* for(i=0; i<=N-1; i+) for(j=0; j<=N-1; j+) A ij=1/(double)(i+j+1); */f0=10; f1=20; f2=30; A00=5; A01=1; A02=2;/ A03=3;A10=1; A11=5; A12=

14、2;/ A13=5;A20=2; A21=1; A22=5;/ A23=3;/ A30=4;/ A31=1;/ A32=2;/ A33=10; for (i=0;i<N;i+) for(j=0;j<N;j+) printf("%f ",A ij); printf("n"); / ldl(A,X,f); printf("n");for(i=0;i<NP;i+) printf("%fn",X i);printf(" LDL's result are:"); for(k=0;

15、k<N;k+) printf("n"); sum=0.0; for(j=0;j<N;j+) sum=sum+Xj*Akj; printf("%f %f", X0,sum); f0=10; f1=20; f2=30; A00=5; A01=1; A02=2;/ A03=3;A10=1; A11=5; A12=2;/ A13=5;A20=2; A21=1; A22=5; gausszy(A,f);/ gauss(A,f); A00=5; A01=1; A02=2;/ A03=3;A10=1; A11=5; A12=2;/ A13=5;A20=2;

16、 A21=1; A22=5; printf("n gauss's result are:"); for(k=0;k<N;k+) printf("n"); sum=0.0; for(j=0;j<N;j+) sum=sum+Xj*Akj; printf("%f, %f", f0,sum); printf("n");解線性方程組之 Gauss消去法#define NP 7#define N 7#include "stdio.h"double XNP;void gauss(doubl

17、e ANPNP,double fNP) int i,j,k; double a_sum=0,temp=0; for(k=0; k<NP; k+) for(i=k+1; i<=NP-1; i+) temp=A ik/Akk; for(j=k+1; j<=NP-1; j+) A ij=A ij-temp*Akj; f i=f i-temp*fk; /回代過程 XNP-1=fNP-1/ANP-1NP-1; for(k=NP-2; k>=0; k-) temp=fk; for(j=NP-1; j>=k+1; j-) temp=temp-Akj*Xj; Xk=temp/Ak

18、k; void main() static double array_aNN; static double array_bN=0,0,0,1,0,0,0; int i,j,k; for(i=0; i<=N-1; i+) for(j=0; j<=N-1; j+) array_a ij=1/(double)(i+j+1); gauss(array_a,array_b); printf(" gauss's result are:"); for(k=0;k<=N-1;k+) printf("n"); printf("%f&qu

19、ot;, Xk); Simpson積分程序設(shè)計(jì)#define N 3 /積分點(diǎn)數(shù)#define a 1.0 /積分下界#define b 2.0 /積分上界#include<iostream.h>int main()float x; float intpN; /積分點(diǎn)坐標(biāo)float wN; /積分系統(tǒng) float f,result; /f表示函數(shù)，result存放積分結(jié)果w0=1.0/6;w1=4.0/6;w2=w0;intp0=a;intp1=(a+b)/2;intp2=b;result=0;for(int i=0;i<N;i+) x=intp i; f=x*x+x+1;

20、result+=w*(b-a)*f;cout<<"積分結(jié)果："<<result<<endl;return 0;牛頓(Newton)插值的算法與程序Newton型插值公式為：Pn (x)=a0+a1(x-x0)+ an(x-x0)(x-x1)(x-xn-1)Newton插值多項(xiàng)式的程序設(shè)計(jì)分兩個(gè)步驟：1是求各階差商an，2是由霍納算法求插值。下面是算法：1. Inputx0,xn, a1,an2. Inputx3. pan4. For i=n-1,n-2,1,0 do 4.1 pp*(x-xi)+ai end do 5. Outputp求a

21、n的算法：1. Input x0,y0,xn,yn2. for i=0 to n do 2.1 ai=yi3. for k=1 to n 3.1. for j=n to k step -1 do k=(ak-ak-1)/(xk-xk-j) end do end do4. Output a0,a1an程序設(shè)計(jì)：#include<iostream.h>float p102; /已知(x0,y0),(x1,y1).，pi 0存xi, pi 1存yifloat a10;int n;int main()int k; float x; /某點(diǎn)插值 float y; /y存放x點(diǎn)對(duì)應(yīng)的插值結(jié)果void an();cout<<"輸入插值組數(shù):"cin>>n;cout<<"輸入"<<n<<"組已知插值數(shù)(X,Y)"<<endl;for(int i=0;i<n;i+) cin>>pi 0>>pi 1;cout<<"輸入插值:"cin>>x;an();y=an-1;for(k=n-2;k&g