- 1、本文档共35页,可阅读全部内容。
- 2、原创力文档(book118)网站文档一经付费(服务费),不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
- 3、本站所有内容均由合作方或网友上传,本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺!文档内容仅供研究参考,付费前请自行鉴别。如您付费,意味着您自己接受本站规则且自行承担风险,本站不退款、不进行额外附加服务;查看《如何避免下载的几个坑》。如果您已付费下载过本站文档,您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“版权申诉”(推荐),也可以打举报电话:400-050-0827(电话支持时间:9:00-18:30)。
查看更多
典型数值算法的C++语言程序设计
1.经典四阶龙格库塔法解一阶微分方程组
1.1算法说明
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。开始,利用下面的计算方法生成近似序列
1.2经典四阶龙格库塔法解一阶微分方程组算法流程图
1.3经典四阶龙格库塔法解一阶微分方程组程序调试
将编写好的代码放在VC6.0环境中编译,直接执行程序便可以得到求解微分方程,并且的结果。如图:
将这些点进行插值或者拟合后就可以得到微分方程的解。
1.4 经典四阶龙格库塔法解一阶微分方程组代码
#include iostream
#include iomanip
using namespace std;
//f为函数的入口地址,x0、y0为初值,xn为所求点,step为计算次数
double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, long step )
{
double k1,k2,k3,k4,result;
double h=(xn-x0)/step;
if(step=0)
return(y0);
if(step==1)
{
k1=f(x0,y0);
k2=f(x0+h/2, y0+h*k1/2);
k3=f(x0+h/2, y0+h*k2/2);
k4=f(x0+h, y0+h*k3);
result=y0+h*(k1+2*k2+2*k3+k4)/6;
}
else
{
double x1,y1;
x1=xn-h;
y1=Runge_Kuta(f, x0, y0, xn-h,step-1);
k1=f(x1,y1);
k2=f(x1+h/2, y1+h*k1/2);
k3=f(x1+h/2, y1+h*k2/2);
k4=f(x1+h, y1+h*k3);
result=y1+h*(k1+2*k2+2*k3+k4)/6;
}
return(result);
}
int main()
{
double f(double x, double y);
double x0,y0;
double a,b;
cout请输入初值x0,y0:;
cinx0y0;
cout请输入区间:;
cinab;
//double x0=0,y0=1;
double x,y,step;
long i;
cout请输入步长:;
cinstep;
//step=0.1;
cout.precision(10);
for(i=0;i=(b-a)/step;i++)
{
x=x0+i*step;
coutsetw(8)xsetw(18)Runge_Kuta(f,x0,y0,x,i)endl;
}
return(0);
}
double f(double x, double y)
{
double r;
r=(x-y)/2;
return(r);
}
2. 高斯列主元法解线性方程组
2.1算法说明
首先将线性方程组做成增光矩阵,对增广矩阵进行行变换。
对第元素,在第i列中,第i行及以下的元素选取绝对值最大的元素,将该元素最大的行与第i行交换,然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。
再对得到的上三角矩阵进行回代操作,即可以得到方程组的解。
2.2高斯列主元算法流程图
2.3高斯列主元程序调试
对所编写的高斯列主元程序进行编译和链接,然后执行得如下所示的窗口,我们按命令输入增广矩阵的行数为4,输入4行5列的增广矩阵:
按回车键后,程序执行得如下所示的结果:
2.4 高斯列主元算法代码
#includeiostream
#includecstdlib
#includecmath
using namespace std;
//在列向量中寻找绝对值最大的项,并返回该项的标号
int FindMax(int p,int N,double *A)
{
int i=0,j=0;
double max=0.0;
for(i=p;iN;i++)
{
if(fabs(A[i*(N+1)+p])max)
{
j=i;
max=fabs(A[i*(N+1)+p]);
}
}
return j;
}
//交换矩阵中的两行
void ExchangeRow(int p,int j,double *A,int N)
{
int i=0;
double C=0
文档评论(0)