1、数值计算课程设计11.经典四阶龙格库塔法解一阶微分方程组1.1 运用四阶龙格库塔法解一阶微分方程组算法分析， ),(1kyxtf)2,2112 ghfhtfkkk,(3 yfxtfkkk),334 hgfhtfkkk),(1kyxtg)2,2112 ghfhtgkkk,(3 yfxtkkk),334 hgfhtgkkk)2(64311 gghyffxk 1kt经过循环计算由 推得 每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为 ,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计NOh算。4 阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精
2、准，稳定，且易于编程。0,txy2,txy（1-1）（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）经典四阶龙格库塔法解一阶微分方程组21.2 经典四阶龙格库塔法解一阶微分方程流程图图 1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3 经典四阶龙格库塔法解一阶微分方程程序代码：#include #include using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double
3、 initial3, double resu3,double h)double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial0;x0=initial1;y0=initial2;f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, 数值计算课程设计3x0+h*f2/2,y0+h*g2/2);f4=f(t
4、0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3);x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu0=t0+h;resu1=x1;resu2=y1;int main()double f(double t,double x, double y);double g(double t,double x, double y);double initial3,resu3;double a,b,H;double t,step;int i;coutinitial0initial1in
5、itial2;coutab;coutstep;cout#include #define N 3using namespace std;void main()int i,j,k,n,p;float t,s,m,aNN,bN,xN;coutaij;coutbi;for(j=0;jfabs(apj)p=i; if(p!=j) /交换第 p 行与第 j 行for(k=0;k=0;i-) s=0.0;for(j=N-1;ji;j-)s=s+aij*xj;xi=(bi-s)/aii; cout#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0
6、.0001 / 差向量 1 范数的上限#define Max 100 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量 yyNvoid ffjacobian(float xxN,float yyNN);/计算雅克比矩阵 yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵 invvoid newdundiedai(float x0N, float invNN,float y0N,float
7、 x1N);/由近似解向量 x0 计算近似解向量 x1float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errorno高斯列主元法解线性方程组12rm;int i,j,iter=0;/如果取消对 x0 的初始化，撤销下面两行的注释符，就可以由键盘向 x0 读入初始近似解向量/for( i=0;ix0i;cout0;i-) for (k=i-1;k=0;k-)L=-augki/augii;for(j=N2-1;j=0;j-)augkj=augkj+L*augij;数值计算课程设计15for (i=0;i=0;i-)for(j=N2-1;j=0
8、;j-)augij=augij/augii;for (i=0;i=k 的近似积分结果逼近表，并以R(j+1,j+1)为最终解来逼近积分。R(j,0)=T(j), j=0,T(j)为区间逐次减半递推梯形求积分公式算出的结果；R(j,1)=S(j), j=1,S(j)为区间逐次减半递推辛普森求积分公式算出的结果；R(j,2)=B(j), j=2,B(j)为递推布尔求积分公式算出的结果；生成 的逼近表 ，并以 为最终解来逼近积分JK(,)RJK(1,)J)(,bafxdR逼近 存在于一个特别的下三角矩阵中，第 0 列元素 用基于 个(,)RJ (,0)RJ2Ja,b子区间的连续梯形方法计算，然后利用
9、龙贝格公式计算 。当K时，第 行的元素为1K(,1)(,1)(,)(,1)4KRJJRJK当 时，程序在第 行结束。|(,)1,|tol()4.2 龙贝格积分算法流程图4)1,(,), kjjkjj （4-1）（4-2）（4-3）数值计算课程设计19图 4-1 龙贝格积分算法流程图4.3 龙贝格积分算法程序代码#includeusing namespace std;#include#include#define f(x) (sin(x) /列举函数龙贝格积分20#define N_H 20#define MAX 10 /最大迭代次数#define a 0 /所求积分的上下限#define b
10、1#define epsilon 0.0001 /所需精度double Romberg(double aa,double bb,long int n)int i;double sum,h=(bb-aa)/n;sum=0;for(i=1;i #include#define MAX 4 double *diff(double X,int n) int i=0; double *H=NULL; H=(double*)malloc(n-1)*sizeof(double); for(i=1;i=1;k-) Mk=(Uk-1-Ck-1*Mk+1)/Bk-1; M0=3*(D0-dx0)/H0-M0/2;
11、MN=3*(dxn-DN-1)/HN-1-MN-1/2; for(k=0;k x1=0:.01:1; y1=polyval(S(1,:),x1-X(1);x2=1:.01:2; y2=polyval(S(2,:),x2-X(2);x3=2:.01:3; y3=polyval(S(2,:),x3-X(3);plot(x1,y1,x2,y2,x3,y3,X,Y,.)运行结果如下：三次样条插值26图 5-2 三次样条插值作图 数值计算课程设计276.M 次多项式曲线拟合6.1 M 次多项式曲线拟合算法说明：设 有 N 个点，横坐标是确定的。最小二乘法抛物线系数表示为kyx1),(，CBxAxfy2)
12、(求解 A，B 和 C 的线性方程组为 43221111322111NNNNkkkkkkkkNNNkkkxxxyxABCxxy整体上考虑近似函数 同所给数据点 (i=0,1,m)误差)(p),(i(i=0,1,m) 的大小，常用的方法有以下三种：一是误差iiiyxpr)(i=0,1,m)绝对值的最大值，即误差向量 的rmr).,(10范数；二是误差绝对值的和，即误差向量 r 的 1范数；三是误差平方和的算术平方根，即误差向量 r 的 2范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2范数的平方，因此在曲线拟合中常采用误差平方和来 度量误差 (i=0，1，m)的整体大小。
13、i数据拟合的具体作法是：对给定数据 (i=0,1,，m)，在取定的函数类),(iyx中,求 ,使误差 (i=0,1,m)的平方和最小。)(xpiiipr)（6-1）（6-2）M 次多项式曲线拟合286.2 M 次多项式曲线拟合流程图FT开始输入 x0、 0Axff)(xfT1)(x?x输出极小值点 和极小值x)(x结束图 6-1 M 次曲线多项式曲线拟合流程图6.3 M 次多项式曲线拟合算法程序代码#include #include #define MAX 20 void inv(double XMAXMAX,int n,double EMAXMAX) int i=0,j=0,k=0; double temp=0.0; 数值计算课程设计29for(i=0;i=0;i-) printf(“%lft“,Ci); printf(“n 拟合后的%d 次多项式为:n“,M); printf(“nP(x)=“); for(i=M;i=0;i-) if(i=0) printf(“%+.3lfn“,Ci); else printf(“%+.3lf*x%d“,Ci,i); getchar(); 6.4 M 次多项式曲线拟合算法程序调试结果：求解（-3，3） ， （0，1） ， （2，1） ， （4，3）四个点拟合 2 次后的 2 后的多项式
