数值课程设计

更新时间：2024-01-01 04:00:01 阅读量：教育文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

数值计算课程设计

1 经典四阶龙格库塔法解一阶微分方程组

1.1 四阶龙格—库塔算法简述

龙格—库塔方法适用于一般的应用，因为它非常精确、稳定，且易于编程。龙格—库塔算法的数学描述如下: 用公式

xk?1?xk?h(f1?2f2?2f3?f4)/6 yk?1?yk?h(g1?2g2?2g3?g4)/6

计算[a,b]上初值问题y'?f(t,y)，y(a)?y0的近似解。

f1?f(tk,xk,yk), g1?g(tk,xk,yk)

hhhhhh,xk?f1,yk?g1) g2?g(tk?,xk?f1,yk?g1) 222222hhhhhhf3?f(tk?,xk?f2,yk?g2) g3?g(tk?,xk?f2,yk?g2)

222222f2?f(tk?f4?f(tk?h,xk?hf3,yk?hg3) g4?g(tk?h,xk?hf3,yk?hg3) 1.2 经典四阶龙格库塔法解一阶微分方程算法流程图

图1-1 四阶龙格库塔算法流程图

四阶龙格—库塔算法

1.3 四阶龙格—库塔算法程序 #include #include using namespace std; #define N 10 int main () {

double fun(double,double,double); double myfun(double,double,double); double f[4],g[4];

double h,a,b;

double x[N+1],y[N+1],t[N+1]; cout<<\输入区间左右端点a,b:\ cin>>a>>b;

cout<

h=(b-a)/N; //h为步长 x[0]=6;

y[0]=4; t[0]=a;

for (int k=0;k

f[1]=fun(t[k]+h/2,x[k]+h/2*f[0],y[k]+h/2*g[0]); g[1]=myfun(t[k]+h/2,x[k]+h/2*f[0],y[k]+h/2*g[0]);

f[2]=fun(t[k]+h/2,x[k]+h/2*f[1],y[k]+h/2*g[1]); g[2]=myfun(t[k]+h/2,x[k]+h/2*f[1],y[k]+h/2*g[1]); f[3]=fun(t[k]+h,x[k]+h*f[2],y[k]+h*g[2]); g[3]=myfun(t[k]+h,x[k]+h*f[2],y[k]+h*g[2]); x[k+1]=x[k]+h*(f[0]+2*f[1]+2*f[2]+f[3])/6; y[k+1]=y[k]+h*(g[0]+2*g[1]+2*g[2]+g[3])/6; t[k+1]=t[k]+h;

}

cout<

数值计算课程设计

return 0; }

double fun(double t,double x,double y) {double f; f=x+2*y; return f; }

double myfun(double t,double x,double y) {double f; f=3*x+2*y; return f; }

1.4 四阶龙格—库塔算法求解例题

在区间[0.0,0.2]上初值x的条件下，求解微分方程组。 6,y40?0?dx?x?2ydtdy?3x?2y dt1.5 四阶龙格—库塔运行结果输入区间左右端点a,b: 0.0 0.2

t x y 0 12.8604 0.22 0.02 6.29355 4.53932 0.04 6.61562 5.11949 0.06 6.96853 5.74397 0.08 7.35474 6.41653 0.1 7.77697 7.14127 0.12 8.23814 7.9226 0.14 8.74141 8.76532 0.16 9.29021 9.6746 0.18 9.88827 10.6561 0.2 10.5396 11.7158

y?11.7158。求得数值解为：x?10.5396 3

高斯列主元法解线性方程组

2 高斯列主元法解线性方程组

2.1 高斯列主元法解线性方程组算法说明

首先将线性方程组做成增广矩阵，对增广矩阵进行行变换。对于元素ai,j，在第i列中，第i行以及以下的元素中选取绝对值最大的元素，将最大元素所在的行与i行交换，然后采取高斯消元法用新得到的ai,j所在的第i行消去第i行以下的元素。依次类推，直到an,n，从而得到上三角矩阵。再进行回代过程。 2.2 高斯列主元法解线性方程组流程图

图2-1 高斯列主元法解线性方程组流程图

数值计算课程设计

2.3 高斯列主元法解线性方程组程序代码 #include #include using namespace std; #define n 4 int main () {

int i,j,k,p,m;

double A[n][n+1],x[n]; double q,L,temp,max1;

cout<<\请输入增广矩阵：\for(i=0;i

for (k=0;k

max1=A[k][k];

for (p=k;p

for(i=0;i

temp=A[k][i]; A[k][i]=A[m][i]; A[m][i]=temp;

for(j=0;j

cin>>A[i][j];

if (fabs(A[p][k])>fabs(max1)) { }

max1=A[p][k];

m=p;

高斯列主元法解线性方程组

}

for(i=k+1;i

L=A[i][k]/A[k][k]; for(j=k;j

{ A[i][j]=A[i][j]-L*A[k][j];

}

cout<<\输出进行上三角迭代过程中的增广矩阵： for(i=0;i

cout<

}

cout<<\输出上三角矩阵：\for(i=0;i cout<

}

x[n-1]=A[n-1][n]/A[n-1][n-1];//回代过程 for(k=n-2;k>=0;k--) { q=0;

for(j=n-1;j>k;j--) { q=q+A[k][j]*x[j];

}

数值计算课程设计

}

x[k]=(A[k][n]-q)/A[k][k];

cout<<\方程组的解为：\for(i=0;i

{

cout<

}

return 0;

2.4 高斯列主元法解线性方程组求解例题

求解下列线性方程组

x1?2x2?72x1?3x2?x3?94x2?2x3?4x4?102x3?4x4?122.5 高斯列主元法解线性方程组算法运行结果请输入增广矩阵： 1 2 0 0 7 2 3 -1 0 9 0 4 2 3 10 0 0 2 -4 12

输出进行上三角迭代过程中的增广矩阵： 2 3 -1 0 9 0 0.5 0.5 0 2.5 0 4 2 3 10 0 0 2 -4 12

输出进行上三角迭代过程中的增广矩阵： 2 3 -1 0 9 0 4 2 3 10

0 0 0.25 -0.375 1.25 0 0 2 -4 12

输出进行上三角迭代过程中的增广矩阵：

高斯列主元法解线性方程组

2 3 -1 0 9 0 4 2 3 10 0 0 2 -4 12 0 0 0 0.125 -0.25 输出上三角矩阵： 2 3 -1 0 9 0 4 2 3 10 0 0 2 -4 12 0 0 0 0.125 -0.25 方程组的解为： 1 3 2 -2

Press any key to continue 求得的解为：x?[1 2 3 ?2]

数值计算课程设计

3 牛顿法解非线性方程组

3.1 牛顿法解非线性方程组算法说明设Pk已知。第一步：计算函数

第二步：计算雅克比矩阵

?f1(pk,qk)?F(Pk)???

f(p,q)?2kk?????f(p,q) f(p,q)1kk1kk??x??y? J(Pk)??????f(p,q) f(p,q)??x2kk?y2kk???第三步：求线性方程组

J(KK)?P??F(PK)

的解?P。

第四步：计算下一点

Pk?1?Pk??P

重复上述过程。

3.2 牛顿法解非线性方程组算法流程图

图3-1 牛顿法解非线性方程组算法流程图 9

牛顿法解非线性方程组

3.3 牛顿法解非线性方程组程序代码 #include #include #include

#define N 2 //用来设置方程组的行数 #define eps 2.2204e-16

double* MatrixMultiply(double* J,double Y[ ]); double* Inv(double *J); double norm(double Q[ ]); double* F(double X[ ]); double* JF(double X[ ]);

int method(double* Y,double epsilon);

int newdim(double P[],double delta,double epsilon,int max1,double *err) {

double *Y=NULL,*J=NULL,Q[2]={0},*Z=NULL,*temp=NULL; double relerr=0.0; int k=0,i=0,iter=0; Y=F(P);

for(k=1;k

J=JF(P);

temp=MatrixMultiply(Inv(J),Y); for(i=0;i

Q[i]=P[i]-temp[i]; Z=F(Q); for(i=0;i

temp[i]=Q[i]-P[i]; *err=norm(temp); relerr=*err/(norm(Q)+eps); for(i=0;i

P[i]=Q[i]; Y[i]=Z[i]; for(i=0;i

数值计算课程设计

}

if((*err

break;

return iter;

int method(double* Y,double epsilon) { }

//矩阵乘法，要求，J为方阵，Y为与J维数相同的列向量 double *MatrixMultiply(double* J,double Y[]) { }

//二阶矩阵的求逆（在M次多项式曲线拟合算法文件中给出了对任意可逆矩阵的求逆算法）

double *Inv(double *J) {

double X[4]={0},temp=0.0; int i=0; X[0]=J[3]; X[1]=-J[1]; double *X=NULL; int i=0,j=0;

X=(double*)malloc(N*sizeof(double)); for(i=0;i

X[i]=0; for(j=0;j

X[i]+=J[i*N+j]*Y[j];

for(i=0;i

if(fabs(Y[0])

return 1; return 0; else

return X;

temp=1/(J[0]*J[3]-J[1]*J[2]);

牛顿法解非线性方程组

}

X[2]=-J[2]; X[3]=J[0]; for(i=0;i<4;i++)

J[i]=temp*X[i]; return J;

double norm(double Q[]) { }

double* F(double X[]) { }

double* JF(double X[]) {

double x=X[0]; double y=X[1]; double *W=NULL;

W=(double*)malloc(4*sizeof(double)); W[0]=2*x-2; double x=X[0]; double y=X[1]; double *Z=NULL;

Z=(double*)malloc(2*sizeof(double)); Z[0]=x*x-2*x-y+0.5; Z[1]=x*x+4*y*y-4; return Z; double max=0.0; int i=0; { }

return max;

if(Q[i]>max)

max=Q[i];

for(i=0;i

数值计算课程设计

}

W[1]=-1; W[2]=2*x; W[3]=8*y; return W;

main() { }

3.4 牛顿法解非线性方程组求解例题

求解非线性方程组

2??x?2x?y?0.5?0 ?22??x?4y?4?0double P[2]={0};

double delta=0.0,epsilon=0.0,err=0.0; int max1=0,iter=0,i=0;

printf(\牛顿法解非线性方程组:\\nx^2-2*x-y+0.5=0\\nx^2+4*y^2-4=0\\n\printf(\输入的初始近似值x0,y0\\n\for(i=0;i<2;i++)

scanf(\

printf(\请依次输入P的误差限，F(P)的误差限，最大迭代次数\\n\scanf(\iter=newdim(P,delta,epsilon,max1,&err); printf(\收敛到P的解为:\\n\for(i=0;i<2;i++)

printf(\printf(\迭代次数为:%d\printf(\误差为:%lf\\n\getchar();

3.5 牛顿法解非线性方程组算法运行结果牛顿法解非线性方程组: x^2-2*x-y+0.5=0 x^2+4*y^2-4=0

牛顿法解非线性方程组

输入的初始近似值x0,y0 2.00 0.25

请依次输入P的误差限，F(P)的误差限，最大迭代次数 0.0001 0.0001 1000 收敛到P的解为: X(1)=1.900691 X(2)=0.311213

迭代次数为:2 误差为:0.000000

20?x?2x?y?0.5非线性方程组在初始值(的条件下求解结p,q)?(2.00,0.25)00220?x?4y?4果为x=1.900691,y=0.311213

数值计算课程设计

4 龙贝格求积分算法

4.1 龙贝格积分算法说明

生成J?K的逼近表R(J,K)，并以R(J?1,J?1)为最终解来逼近积分逼近R(J,K)存在一个特别的三角矩阵中，第0列元素R(J,0)用基于2J个?a,b?子区间的连续梯形方法计算，然后利用龙贝格公式计算R(J,K)。当1?K?J时，第J行的元素为

R(J,K?1)?R(J?1,K?1)

4K?1当R(J,J)?R(J?1,J?1)?tol时，程序在第?J?1?行结束。

R(J,K)?R(J,K?1)??baf(x)dx?R(J,J)

4.2 龙贝格求积分算法的算法流程图

图4-1 龙贝格求积分的算法流程图

4.3 龙贝格积分算法程序代码

龙贝格积分算法

#include #include #include using namespace std; #define n 8 int main() {

double feval(double);//函数声明 double a=0,b=3,h,x,s,err=1; double R[n][n]={0}; int M=1,J=0; h=b-a;

R[0][0]=h*(feval(a)+feval(b))/2; while ((err>0.01)&&(J

cout<<\该矩阵为:\for (int i=0;i

J=J+1; h=h/2; s=0;

for (int p=1;p

for (int K=1;K<=J;K++) { }

err=fabs(R[J-1][J-1]-R[J][K]);

R[J][K]=(pow(4.0,K)*R[J][K-1]-R[J-1][K-1])/(pow(4.0,K)-1); x=a+h*(2*p-1); s=s+feval(x);

R[J][0]=R[J-1][0]/2+h*s;

数值计算课程设计

}

for (int j=0;j<=i;j++)

cout<

cout<<\则积分为=\ cout<

double feval(double x) //定义函数 { double f;

f=sin(2*x)/(1+x*x); return f;

}

4.4 龙贝格积分求解例题

求?3sin(2x)91?x2的值。 4.5 龙贝格积分算法运行结果该矩阵为:

-0.041912325 0.044176149 0.072872307 0.37995411 0.49188009 0.51981395 0.45457605 0.47945003 0.47862136 0.47796751 0.47086479 0.47629437 0.47608399 0.47604371 0.47603617 0.4748323 0.47615481 0.4761455 0.47614648 0.47614688 0.47581819 0.47614682 0.47614629 0.4761463 0.4761463 0.4760643 0.47614633 0.4761463 0.4761463 0.4761463 则积分为= 0.4761463 Press any key to continue

0.47614699 0.4761463 0.4761463 0.4761463 0.4761463 0.4761463 三次样条插值算法

5 三次样条插值算法（压紧样条）

5.1 三次样条插值算法（压紧样条）说明

表5-1 三次样条插值算法（压紧样条）说明

策略描述 (i)三次紧压样条，确定S?(x0)，S?(xn)（如果导数已知，这是“最佳选择”）包含m0和mN的方程 m0?m3(d0?S?(x0))?1 h023hN?1(S?(xN)?dN?1)?mN?1 2mN?(ii)natural三次样条（一条“松弛曲线”） (iii)外挂S??(x)到端点 m0?0，mN?0 m0?m1?h0(m2?m1) h1mN?mN?1?S??(x)是靠近端点的常量 hN?1(mN?1?mN?2) hN?2 (iv) m0?m1，mN?mN?1(v)在每个端点处指定S??(x) m0?S??(x0)，mN?S??(xN) 5.2 三次样条插值算法（压紧样条）的算法流程图

开始计算S[k][0]S[i][1],S[i][2],S[i][3] 输出S[i][0],S[i][1],S[i][2],S[i][3] 结束图5-1 三次样条插值算法（压紧样条）的算法流程 18

数值计算课程设计

5.3 三次样条插值算法（压紧样条）程序代码 #include #include #define MAX 4

double *diff(double X[],int n) { }

double *divide(double Y[],int N,double H[]) { } main() {

double

X[MAX]={0,1,2,3},Y[MAX]={0,0.5,2.0,1.5},S[MAX][MAX]={0},temp=0.0,M[MAX]={0};

int i=0;

double *H=NULL;

H=(double*)malloc((n-1)*sizeof(double)); for(i=1;i<=n-1;i++) { } return H;

H[i-1]=X[i]-X[i-1];

int i=0;

double *D=NULL;

D=(double*)malloc(N*sizeof(double)); for(i=0;i

D[i]=Y[i]/H[i];

int N=MAX-1,i=0,k=0;

double A[MAX-1-2]={0},B[MAX-1-1]={0},C[MAX-1-1]={0}; double dx0=0.2,dxn=1.0;

三次样条插值算法

double *H=NULL,*D=NULL,*U=NULL;

printf(\求解经过点(0,0.0),(1,0.5),(2,2.0)和(3,1.5)\\n而且一阶导数边界条件

S'(0)=0.2和S'(3)=-1的三次压紧样条曲线\\n\\n\H=diff(X,MAX);

D=divide(diff(Y,MAX),N,H); for(i=1;i

A[i]=H[i+1]; B[i]=2*(H[i]+H[i+1]); C[i]=H[i+1]; for(i=0;i

U[i]=U[i]*6; B[0]=B[0]-H[0]/2; U[0]=U[0]-3*(D[0]-dx0); B[N-2]=B[N-2]-H[N-1]/2; U[N-2]=U[N-2]-3*(dxn-D[N-1]); for(k=2;k<=N-1;k++) {

temp=A[k-2]/B[k-2]; B[k-1]=B[k-1]-temp*C[k-2]; U[k-1]=U[k-1]-temp*U[k-2]; }

M[N-1]=U[N-2]/B[N-2]; for(k=N-2;k>=1;k--)

M[k]=(U[k-1]-C[k-1]*M[k+1])/B[k-1]; M[0]=3*(D[0]-dx0)/H[0]-M[0]/2; M[N]=3*(dxn-D[N-1])/H[N-1]-M[N-1]/2; {

S[k][0]=(M[k+1]-M[k])/(6*H[k]); S[k][1]=M[k]/2;

S[k][2]=D[k]-H[k]*(2*M[k]+M[k+1])/6;

for(k=0;k<=N-1;k++)

数值计算课程设计

}

S[k][3]=Y[k];

cout<<\输出s[N][4]:\

for(i=0;i

cout<

for(k=0;k<4;k++)

cout<

5.4 三次样条插值算法（压紧样条）求解例题

求三次紧压样条曲线，我们以经过点（0,0），（1,0.5），(2,2.0)和（3,1.5），

‘‘且一阶导数的边界条件为S(0)?0.2和S(3)?-1的拟合函数

5.5 三次样条插值算法（压紧样条）程序运行结果

输出s[N][4]:

0.48 -0.18 0.2 0 -1.04 1.26 1.28 0.5 0.68 -1.86 0.68 2

5.6 运用Matlab绘制三次压紧样条的函数图像我们借助Matlab绘制出以上三次压紧样条的函数图像

图5-2 三次压紧样条的函数图像

M次多项式曲线拟合

6 M次多项式曲线拟合

6.1 M次多项式曲线拟合算法说明

设?有N个点，横坐标是确定的。最小二乘抛物线的系数表示为 (x?k,yk)k?12y?fx()?Ax?BxC?

N求解A，B和C的线性方程组为

N?N4??N3??N2?2??xk?A???xk?B???xk?C??ykxkk?1?k?1??k?1??k?1?

N?N3??N2??N???xk?A???xk?B???xk?C??ykxkk?1?k?1??k?1??k?1?N?N2??N?C??yk??xk?A???xk?B?Nk?1?k?1??k?1?

6.2 M次多项式曲线拟合算法流程图

图6-1 算法流程图

6.3 M次多项式曲线拟合算法程序 #include

数值计算课程设计

#include #include #define n 6 #define M 4 using namespace std; int main() {

int i,j,k,p,m;

double X[1][n],Y[1][n],F[n][M+1];

double A[M+1][M+1]={0},B[M+1][1]={0},C[M+1][1]; double E[M+1][M+1],D[M+1][2*M+2],AN[M+1][M+1];

//E[M+1][M+1]为单位矩阵,D[M+1][2*M+2]矩阵形式为(A,E),AN[M+1][M+1]为A[M+1][M+1]的逆矩阵。

double XZ[n][1],YZ[n][1],FZ[M+1][n];

// XZ[n][1]为X[1][n]矩阵的转置，YZ[n][1]为Y[n][1]矩阵的转置，FZ[M+1][n]为F[n][M+1]矩阵的转置. double temp,L,max1;

cout<<\输入矩阵X[1][n]:\for (i=0;i

cin>>X[0][i];

cout<<\输入矩阵Y[1][n]:\for (i=0;i

cin>>Y[0][i];

cout<<\输出矩阵XZ[1][n]:\求矩阵XZ[1][n] for (i=0;i

cout<<\输出矩阵YZ[1][n]:\求矩阵YZ[1][n] for (i=0;i

YZ[i][0]=Y[0][i]; XZ[i][0]=X[0][i];

cout<

M次多项式曲线拟合

}

cout<

cout<<\输出矩阵F[n][M+1]:\计算矩阵F[n][M+1] for(k=0;k

for(i=0;i

cout<<\输出矩阵FZ[M+1][n]:\计算矩阵FZ[M+1][n] for(i=0;i

FZ[k][i]=F[i][k]; }

for(k=0;k

for (i=0;i

for (j=0;j

cout<

F[i][k]=pow(XZ[i][0],k);

for(k=0;k

cout<

for(k=0;k

for(i=0;i

cout<

cout<<\输出矩阵A[M+1][M+1]=\计算A[M+1][M+1]

数值计算课程设计

{

}

for (k=0;k

A[i][j]=A[i][j]+FZ[i][k]*F[k][j];

for (i=0;i

for (j=0;j

cout<}

for (i=0;i

cout<

B[i][0]=B[i][0]+FZ[i][j]*YZ[j][0];

cout<<\输出矩阵B[M+1][1]=\计算B[M+1][1]

cout<

cout<<\输出矩阵E[M+1][M+1]=\构造单位矩阵 for (i=0;i

for (i=0;i

for (j=0;j

if (i==j)

E[i][j]=1; E[i][j]=0; else

M次多项式曲线拟合

{

for (j=0;j

cout<

}

cout<<\输出矩阵D[M+1][2*M+2]=\构造形式为(A,E)的矩阵for (i=0;i

{

for (j=0;j<2*M+2;j++)

cout<

for (k=0;k

max1=D[k][k];

for (p=k;p

for (i=0;i<2*M+2;i++)//交换行 {

temp=D[k][i]; D[k][i]=D[m][i];

if (fabs(D[p][k])>fabs(max1)) { }

max1=D[p][k]; m=p;

for (j=0;j

D[M+1][2*M+2]

for (i=0;i

for (j=M+1;j<2*M+2;j++)

D[i][j]=E[i][j-M-1];

for (i=0;i

数值计算课程设计

D[m][i]=temp; {

}

for(i=k-1;i>=0;i--) {

L=D[i][k]/D[k][k]; for(j=2*M+1;j>=0;j--)

D[i][j]=D[i][j]-L*D[k][j];

} {

for (i=k+1;i

L=D[i][k]/D[k][k]; }

for(k=M;k>0;k--)

D[i][j]=D[i][j]-L*D[k][j];

for (j=k;j<2*M+2;j++) }

}

for (i=M;i>=0;i--) //将矩阵D[M+1][2*M+2]中的A[M+1][M+1]变为单位矩阵

{

for (j=2*M+1;j>=0;j--)

D[i][j]=D[i][j]/D[i][i];

}

cout<<\输出最终的D[M+1][2*M+2]:\

cout<<\输出矩阵AN[M+1][M+1]=\计算矩阵AN[M+1][M+1] for (i=0;i

{

for (j=0;j

AN[i][j]=D[i][j+M+1];

}

for (i=0;i

for (j=0;j

M次多项式曲线拟合

}

cout<cout<

cout<<\输出矩阵C[M+1][1](矩阵C[M+1][1]=A/B,即C=AN*B)=\ for (i=0;i

6.4 M次多项式曲线拟合算法求解例题

通过点（0.25,23.1），（1.0,1.68）,(1.5,1.0),(2.0,0.84),(2.4,0.826),(5.0,1.2576)拟合一个四次多项式。

6.5 M次多项式曲线拟合算法运行结果输入矩阵X[1][n]: 0.25 1.0 1.5 2.0 2.4 5.0 输入矩阵Y[1][n]:

23.1 1.68 1.0 0.84 0.826 1.2576 输出矩阵XZ[1][n]:

{ }

cout<<\输出矩阵(对矩阵C[M+1][1])进行翻转C[M+1][1]=\{ }

for (i=0;i

cout<

temp=C[M-i][0]; C[i][0]=temp; C[i][0]=0;

for (j=0;j

C[i][0]=C[i][0]+AN[i][j]*B[j][0]; cout<

for (i=0;i<(M+1)/2;i++)

C[M-i][0]=C[i][0];

数值计算课程设计

0.25 1 1.5 2 2.4 5

输出矩阵YZ[1][n]: 23.1 1.68 1 0.84 0.826 1.2576

输出矩阵F[n][M+1]:

1 0.25 0.0625 0.015625 0.00390625 1 1 1 1 1 1 1.5 2.25 3.375 5.0625 1 2 4 8 16 1 2.4 5.76 13.824 33.1776 1 5 25 125 625 输出矩阵FZ[M+1][n]:

1 1 1 1 1 1 0.25 1 1.5 2 2.4 5 0.0625 1 2.25 4 5.76 25 0.015625 1 3.375 8 13.824 125 0.00390625 1 5.0625 16 33.1776 625 输出矩阵A[M+1][M+1]=

6 12.15 38.0725 151.215 680.244 12.15 38.0725 151.215 680.244 3245.22 38.0725 151.215 680.244 3245.22 15892.5 151.215 680.244 3245.22 15892.5 78729.7 680.244 3245.22 15892.5 78729.7 392008 输出矩阵B[M+1][1]=

28.7036 18.9054 44.9315 180.755 833.677

M次多项式曲线拟合

输出矩阵E[M+1][M+1]=

1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 输出矩阵D[M+1][2*M+2]=

6 12.15 38.0725 151.215 680.244 1 0 0 0 0

12.15 38.0725 151.215 680.244 3245.22 0 1 0 0 0

38.0725 151.215 680.244 3245.22 15892.5 0 0 1 0 0

151.215 680.244 3245.22 15892.5 78729.7 0 0 0 1 0

680.244 3245.22 15892.5 78729.7 392008 0 0 0 0 1

输出矩阵AN[M+1][M+1]=

7.60557 -30.2544 31.0739 -11.273 1.24151 -30.2544 139.476 -150.341 55.7933 -6.21249 31.0739 -150.341 166.101 -62.5425 7.01758 -11.273 55.7933 -62.5425 23.7722 -2.68111 1.24151 -6.21249 7.01758 -2.68111 0.303242 输出矩阵C[M+1][1](矩阵C[M+1][1]=A/B,即C=AN*B)= 39.9168 -80.9324 58.3927 -17.1522 1.68035

输出矩阵(对矩阵C[M+1][1])进行翻转C[M+1][1]= 1.68035 -17.1522 58.3927 -80.9324 39.9168 Press any key to continue

得到的拟合函数为：y?1.68035x4?17.1522x3?58.3927x2?80.9324x?39.9168 6.5 Matlab画图并观察拟合效果

数值计算课程设计

x=linspace(-10,400,10);

y=1.68035*x.^4-17.1522*x.^3+58.3927*x.^2-80.9324*x+39.9168; plot(x,y,'b')

4.543.532.521.510.50-50x 1010050100150200250300350400

图6-2 多项式拟合的函数图像

试值法或试位法求解方程的根

7 试值法或试位法求解方程的根

7.1 试值法或试位法算法说明

求解方程f(x)?0在区间[a,b]内的根，前提条件是与f(b)的符号相反。 7.2 试值法或试位法算法程序 #include #include using namespace std; int main ( ) {

且f(a)f(x是连续的，

double feval(double); double min1(double,double); int k,max1=100;

double a,b,c,ya,yb,yc,dx,ac,err,delta=0.0001,epsilon=0.0001; cout<<\输入区间端点值a,b:\cin>>a>>b; ya=feval(a); yb=feval(b);

cout<<\输出方程的根:\for(k=1;k<=max1;k++) {

dx=yb*(b-a)/(yb-ya); c=b-dx; ac=c-a; yc=feval(c); if (yc==0) { }

else if (yb*yc>0) {

b=c; yb=yc;

cout<

数值计算课程设计

}

} else { }

dx=min1(fabs(dx),ac);

if ((fabs(dx)

break; a=c; ya=ya;

err=fabs(b-a)/2; yc=feval(c); cout<

cout<<\cout<<\return 0;

double feval(double x) { }

double min1(double dx,double ac) { }

if (fabs(dx)>ac)

dx=ac; else dx=fabs(dx); return 0;

return(x*sin(x)-1);

7.3 试值法或试位法算法求解例题

求解xsin(x)?1?0在区间[0,2]内的根。 7.4 试值法或试位法算法运行结果输入区间端点值a,b:0 2 输出方程的根:1.09975

试值法或试位法求解方程的根

err=:0.450125 yc=:-0.0200192 Press any key to continue

数值计算课程设计

8 牛顿—拉夫森迭代算法

8.1 牛顿—拉夫森迭代算法说明使用初始近似值p0利用迭代

pk?pk?1?f(pk?1) 其中k=1,2,3,......

'f(pk?1)计算函数f(x)?0的根的近似值。 8.2 牛顿—拉夫森迭代算法程序代码 #include #include using namespace std; int main () { double f1(double); double df1(double); int k,max1=100;

double p0,p1,y,err,relerr,delta=0.1,epsilon=0.1; cout<<\ cin>>p0;

cout<<\输出迭代次数:\输出方程的根:\ for (k=1;k<=max1;k++) { p1=p0-f1(p0)/df1(p0); err=fabs(p1-p0);

relerr=2*err/(fabs(p1)+delta); p0=p1; y=f1(p0);

if ((err

cout<

cout<<\输出方程的根:\

牛顿—拉夫森迭代算法

cout<<\输出迭代次数:\ cout<<\ cout<<\ return 0;

}

double f1(double x) { double f; f=pow(x,2.0)-2.0;

return f;

}

double df1(double x) { double df; df=2*x; return df;

}

8.3 牛顿—拉夫森迭代算法求解例题求解x2?2?0的根。

8.4 牛顿—拉夫森迭代算法运行结果p0=:2

输出迭代次数，输出方程的根: 1 1.5 0.25 输出方程的根:1.41667 输出迭代次数:2 err:0.0833333 relerr:0.10989

Press any key to continue

数值计算课程设计

9 利用外推法的微分求解

9.1 外推算法说明

为了求解f'(x)的数值解，构造包含近似值D(j，k)，k?j的表，并将

f'(x)?D（n，n作为最终答案。近似值D(j，k)存放在下三角矩阵中。

第一列是

f(x?2?jh)?f(x?2?jh) D(j,0)??j?12h第j行元素为

D(j,k)?D(j,k?1)?D(j,k?1)?D(j?1,k?1)，其中1?k?j k4?19.2 外推算法程序代码 #include #include #include #define n 4

using namespace std; int main() {

double feval(double); int i,j=1,k,p;

double x,delta=0.0001,toler=0.0001,eps=0.01; double h,err=1,relerr=1,D[n][n]; cout<<\输入x值:\cin>>x;

cout<<\输入h值:\cin>>h;

D[0][0]=(feval(x+h)-feval(x-h))/(2*h); while ((relerr>toler)&&(err>delta)&&(j<12)) {

h=h/2;

D[j][0]=(feval(x+h)-feval(x-h))/(2*h); for (k=1;k<=j;k++) {

利用外推法的微分求解

}

D[j][k]=D[j][k-1]+(D[j][k-1]-D[j-1][k-1])/(pow(4.0,k)-1);

}

err=fabs(D[j][j]-D[j-1][j-1]);

relerr=2*err/(fabs(D[j][j])+fabs(D[j-1][j-1])+eps); }

for ( i=0;i

cout<<\则微分为=\cout<

cout<<\cout<<\return 0;

for (int p=0;p<=i;p++)

cout<

j=j+1;

cout<<\该矩阵为:\

double feval(double x) { }

9.3 外推算法求解例题

求解f(0.8)?sin(0.8)的值。 9.4 外推算法运行结果输入x值:0.8 输入h值:0.01 该矩阵为: -0.71734414

-0.7173531 -0.71735609

double f; f=cos(x); return f;

数值计算课程设计

则微分为= -0.71735609 err=:1.195586e-005 relerr=:1.6551337e-005 Press any key to continue

组合辛普森公式

10 组合辛普森公式

10.1 组合辛普森公式算法说明

通过f(x)的2M+1个等步长采样点xk?a?kh,k?0,1,2?，2M逼近积分

h2hM?14hM?f(x)dx?(f(a)?f(b))?f(x2k?1) ?f(x2k)?3k??133k?1ba注意， x0?a，x2M?b。

10.2 组合辛普森公式程序代码#include #include using namespace std; int main() { double feval(double); int k,M;

double a,b,h,s1,s2,s,x; s1=0; s2=0;

cout<<\输入M值:\ cin>>M;

cout<<\输入左右端点值:\ cin>>a>>b;

h=(b-a)/(2*M); for (k=1;k<=M;k++) { x=a+h*(2*k-1); s1=s1+feval(x);

}

for (k=1;k<=M-1;k++) { x=a+h*2*k; s2=s2+feval(x);

}

数值计算课程设计

}

s=h*(feval(a)+feval(b)+4*s1+2*s2)/3; cout<<\输出s:\return 0;

double feval(double x) { }

10.3 组合辛普森公式求解例题

对函数f(x)?x3在区间[0,1]上积分近似值。 10.4 组合辛普森公式运行结果

double f; f=pow(x,3); return f;

图 10.1 运行结果

设计体会及意见

11. 设计体会及意见

11.1 设计体会

这次课程设计主要是将Maltab语言翻译成C++语言，但是还没有熟练地掌握数值算法以及C++的编程能力还不好，所以不能使编程达到很完美，另外，这次课程设计让我们认识到自己所学的知识有限，技术尚未达到很高的水准，不得不虚心向老师及同学求教，以及多多地查询有关的资料，就这方面来说，它让我明白了自身的不足之处，让我懂得在学习的过程更加虚心！。但是通过这次课程设计的实践既提高了我对Maltab程序的应用，又提高了对C++程序的编程，还提高了对Word的编辑，培养了我的独立思考、分析以及解决问题能力，同样也磨练了我的毅力。 12.1 改进意见

在此次课程设计中，由于大部分是参考Maltab程序，没有按照自己的想法去完成，所以应该多独立思考，首先应熟练的掌握算法，然后在Maltab实现，最后完成设计。此次设计中，在求矩阵的逆时，比较繁琐，需要进一步改进。

参考文献

[1] John H.Matews.数值方法（第四版）[M].北京：电子工业出版社.2010.4 [2] 谭浩强.C++程序设计[M].北京：清华大学出版社.2004.6

1 经典四阶龙格库塔法解一阶微分方程组................................ 1

1.1 四阶龙格—库塔算法简述 ...................................... 1 1.2 经典四阶龙格库塔法解一阶微分方程算法流程图 .................. 1 1.3 四阶龙格—库塔算法程序 ...................................... 2 1.4 四阶龙格—库塔算法求解例题 .................................. 3 1.5 四阶龙格—库塔运行结果 ...................................... 3 2 高斯列主元法解线性方程组.......................................... 4

2.1 高斯列主元法解线性方程组算法说明 ............................ 4 2.2 高斯列主元法解线性方程组流程图 .............................. 4 2.3 高斯列主元法解线性方程组程序代码 ............................ 5 2.4 高斯列主元法解线性方程组求解例题 ............................ 7 2.5 高斯列主元法解线性方程组算法运行结果 ........................ 7 3 牛顿法解非线性方程组.............................................. 9

3.1 牛顿法解非线性方程组算法说明 ................................ 9 3.2 牛顿法解非线性方程组算法流程图 .............................. 9 3.3 牛顿法解非线性方程组程序代码 ............................... 10 3.4 牛顿法解非线性方程组求解例题 ............................... 13 3.5 牛顿法解非线性方程组算法运行结果 ........................... 13 4 龙贝格求积分算法................................................. 15

4.1 龙贝格积分算法说明 ......................................... 15 4.2 龙贝格求积分算法的算法流程图 ............................... 15 4.3 龙贝格积分算法程序代码 ..................................... 15 4.4 龙贝格积分求解例题 ......................................... 17 4.5 龙贝格积分算法运行结果 ..................................... 17 5 三次样条插值算法（压紧样条）..................................... 18

5.1 三次样条插值算法（压紧样条）说明 ........................... 18 5.2 三次样条插值算法（压紧样条）的算法流程图 ................... 18 5.3 三次样条插值算法（压紧样条）程序代码 ....................... 19 5.4 三次样条插值算法（压紧样条）求解例题 ....................... 21 5.5 三次样条插值算法（压紧样条）程序运行结果 ................... 21 5.6 运用Matlab绘制三次压紧样条的函数图像 ...................... 21 6 M次多项式曲线拟合 ............................................... 22

6.1 M次多项式曲线拟合算法说明.................................. 22

6.2 M次多项式曲线拟合算法流程图................................ 22 6.3 M次多项式曲线拟合算法程序.................................. 22 6.4 M次多项式曲线拟合算法求解例题.............................. 28 6.5 M次多项式曲线拟合算法运行结果.............................. 28 6.5 Matlab画图并观察拟合效果................................... 30 7 试值法或试位法求解方程的根....................................... 32

7.1 试值法或试位法算法说明 ..................................... 32 7.2 试值法或试位法算法程序 ..................................... 32 7.3 试值法或试位法算法求解例题 ................................. 33 7.4 试值法或试位法算法运行结果 ................................. 33 8 牛顿—拉夫森迭代算法............................................. 35

8.1 牛顿—拉夫森迭代算法说明 ................................... 35 8.2 牛顿—拉夫森迭代算法程序代码 ............................... 35 8.3 牛顿—拉夫森迭代算法求解例题 ............................... 36 8.4 牛顿—拉夫森迭代算法运行结果 ............................... 36 9 利用外推法的微分求解............................................. 37

9.1 外推算法说明 ............................................... 37 9.2 外推算法程序代码 ........................................... 37 9.3 外推算法求解例题 ........................................... 38 9.4 外推算法运行结果 ........................................... 38 10 组合辛普森公式.................................................. 40

10.1 组合辛普森公式算法说明 .................................... 40 10.2 组合辛普森公式程序代码 .................................... 40 10.3 组合辛普森公式求解例题 .................................... 41 10.4 组合辛普森公式运行结果 .................................... 41 11. 设计体会及意见................................................. 42

11.1 设计体会 .................................................. 42 12.1 改进意见 .................................................. 42 TOC \\o \参考文献....................................... 43

本文来源：https://www.bwwdw.com/article/fpix.html

相关文章：

正在阅读：

数值课程设计01-01

杂多县职称论文发表网-建筑工程消防自动控制系统设计论文选题题04-10

做一个感恩的员工02-18

第四讲 Microsoft PowerPoint 演示文稿08-26

大学英语1复习资料04-05

家庭社会工作考试试题及答案解析12-27

民用建筑工程的绿色施工技术探讨_004-22

ISO 90012015新版标准正式发布05-23

最新-浙江省2018各市地科学中考分册汇编精品12-19

中高级爆破工程技术人员安全作业证04-08

上一篇：浅析美国量化宽松货币政策的影响下一篇：圣训经评论