《[数学]数值计算课程设计》由会员分享，可在线阅读，更多相关《[数学]数值计算课程设计（35页珍藏版）》请在金锄头文库上搜索。
1、数值计算课程设计1、 经典四阶龙格库塔法解一阶微分方程组1.1 经典四阶龙格库塔法解一阶微分方程组算法说明4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点开始，利用生成的近似序列，其中1.2 经典四阶龙格库塔法解一阶微分方程组算法流程图 图1-1第 1 页 共 38 页数值计算课程设计1.3 经典四阶龙格库法解一阶微分方程组程序调试 图1-21.4 经典四阶龙格库法解一阶微分方程组程序代码#include #includeusing namespace std;int main() double f1(double,double); double f2(d
2、ouble,double); int t,j,N=10; double a=0,b=0.2,T12,Z123,h; h=(b-a)/N; for(t=1;t12;t+) Tt=a+(t-1)*h; Z11=6;Z12=4; double k1,k2,k3,k4,m1,m2,m3,m4; for(j=1;jN+1;j+) k1=h*f1(Zj1,Zj2); m1=h*f2(Zj1,Zj2); k2=h*f1(Zj1+k1/2,Zj2+m1/2); m2=h*f2(Zj1+k1/2,Zj2+m1/2); k3=h*f1(Zj1+k2/2,Zj2+m2/2); m3=h*f2(Zj1+k2/2,Zj
3、2+m2/2); k4=h*f1(Zj1+k3,Zj2+m3); m4=h*f2(Zj1+k3,Zj2+m3); Zj+11=Zj1+(k1+2*k2+2*k3+k4)/6; Zj+12=Zj2+(m1+2*m2+2*m3+m4)/6; cout -endl;cout 计算微分方程组： endl;cout dx/dt=3x+5y endl;cout dy/dt=2x+y endl;cout 的龙格-库塔解 endl;cout -endl;cout.setf(ios:fixed); coutsetprecision(5); coutsetw(6)tksetw(15)xksetw(15)ykend
4、l; for(t=1;t12;t+) coutsetw(8)Ttsetw(15)Zt1setw(15)Zt2endl; return 0;double f1(double x,double y)double z;z=x+2*y;return z;double f2(double x1,double y1)double z1;z1=3*x1+2*y1;return z1;2、 高斯列主元法解线性方程组21 高斯列主元法解线性方程组算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元
5、法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2高斯列主元法解线性方程组流程图 图2-12.3高斯列主元法解线性方程组程序调试 图2-22.4高斯列主元法解线性方程组程序代码#include#include #define N 3using namespace std;void main()int i,j,k,n,p;float t,s,m,aNN,bN,xN;cout请输入方程组的系数endl;for(i=0;iN;i+)for(j=0;jaij;cout请输入方程组右端的常数项:endl;for(i=0;ib
6、i;for(j=0;jN-1;j+) p=j; for(i=j+1;ifabs(apj) p=i; if(p!=j) /交换第p行与第j行for(k=0;kN;k+)t=ajk;ajk=apk;apk=t; /交换常数项的第p行与第j行 t=bp; bp=bj; bj=t; /把系数矩阵第j列ajj下面的元素变为0for(i=j+1;iN;i+) m=-aij/ajj; for(n=j;n=0;i-) s=0.0;for(j=N-1;ji;j-)s=s+aij*xj; xi=(bi-s)/aii; cout方程组的解如下:endl;for(i=0;i=N-1;i+) coutxiendl;第
7、35 页 共 37 页数值计算课程设计3、 牛顿法解非线性方程组31 牛顿法解非线性方程组的基本算法设已知。第1步：计算函数 第2步：计算雅可比矩阵第3步：求线性方程组的解。第4步：计算下一点重复上述过程。3.2牛顿法解非线性方程组的流程图 图3-13.3牛顿法解非线性方程组的程序调试 图3-23.4牛顿法解非线性方程组的程序代码#include#include#define N 2#define Epsilon 0.00001#define Max 200 using namespace std;const int N2=2*N;int main()void ff(float xxN,flo
8、at yyN);void ffjacobian(float xxN,float yyNN);void inv_jacobian(float yyNN,float invNN);void newdundiedai(float x0N, float invNN,float y0N,float x1N);float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;cout初始近似解向量：endl;for (i=0;iN;i+) coutx0i ; coutendl;coutendl; doiter=iter
9、+1;cout第 iter 次迭代开始endl;ff(x0,y0); ffjacobian(x0,jacobian); inv_jacobian(jacobian,invjacobian); newdundiedai(x0, invjacobian,y0,x1); errornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break; for (i=0;iN;i+)x0i=x1i; while (iterMax);return 0;void ff(float xxN,float yyN)
10、float x,y; int i; x=xx0;y=xx1; yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4; cout向量函数的因变量向量是： endl;for( i=0;iN;i+) coutyyi ; coutendl;coutendl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j; x=xx0; y=xx1; yy00=2*x-2; yy01=-1; yy10=2*x; yy11=8*y; cout雅克比矩阵是： endl;for( i=0;iN;i+) for(j=0;jN;j+) coutyyij ; coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; cout开始计算雅克比矩阵的逆矩阵 ：endl; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else
