经典四阶龙格库塔法解一阶微分方程组Word文档格式.docx
《经典四阶龙格库塔法解一阶微分方程组Word文档格式.docx》由会员分享,可在线阅读,更多相关《经典四阶龙格库塔法解一阶微分方程组Word文档格式.docx(38页珍藏版)》请在冰豆网上搜索。
y1=y0+h*(g1+2*g2+2*g3+g4)/6;
resu[0]=t0+h;
resu[1]=x1;
resu[2]=y1;
}
intmain()
doublef(doublet,doublex,doubley);
doubleg(doublet,doublex,doubley);
doubleinitial[3],resu[3];
doublea,b,H;
doublet,step;
inti;
cout<
<
"
输入所求微分方程组的初值t0,x0,y0:
;
cin>
>
initial[0]>
initial[1]>
initial[2];
输入所求微分方程组的微分区间[a,b]:
a>
b;
输入所求微分方程组所分解子区间的个数step:
step;
setiosflags(ios:
:
right)<
fixed)<
setprecision(10);
H=(b-a)/step;
initial[0]<
setw(18)<
initial[1]<
initial[2]<
endl;
for(i=0;
i<
i++)
{RK4(f,g,initial,resu,H);
cout<
resu[0]<
setw(20)<
resu[1]<
resu[2]<
initial[0]=resu[0];
initial[1]=resu[1];
initial[2]=resu[2];
return(0);
doublef(doublet,doublex,doubley)
doubledx;
dx=x+2*y;
return(dx);
doubleg(doublet,doublex,doubley)
doubledy;
dy=3*x+2*y;
return(dy);
1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示:
应用所编写程序计算所给例题:
其中初值为
求解区间为[0,0.2]。
计算结果为:
图1-2经典四阶龙格库塔法解一阶微分方程算法程序调试图
2.高斯列主元法解线性方程组
2.1高斯列主元法解线性方程组算法分析
使用伪代码编写高斯消元过程:
fork=1ton-1do
fori=k+1ton
l<
=a(i,k)/a(k,k)
forj=ktondo
a(i,j)<
=a(i,j)-l*a(k,j)
end
%endofforj
b(i)<
=b(i)-l*b(k)
%endoffori
%endoffork
最后得到A,b可以构成上三角线性方程组
接着使用回代法求解上三角线性方程组
因为高斯消元要求a(k,k)≠0(k=1,2,3……n-1)这就需要对高斯消元过程进行完善,即使用高斯列主元法:
其步骤为:
①找主元:
计算
并记录其所在行r,
②交换第r行与第k行;
③以第k行为工具行处理以下各行,使得从第k列的第k+1行到第n行的元素全部为0;
④得到增广矩阵的上三角线性方程组;
⑤使用回代法对上三角线性方程组进行求解
2.2高斯列主元法解线性方程组流程图
图2-1高斯列主元法解线性方程组流程图
2.3高斯列主元法解线性方程组程序代码
#include<
cmath>
#defineN3
voidmain()
{inti,j,k,n,p;
floatt,s,m,a[N][N],b[N],x[N];
请输入方程组的系数"
N;
{for(j=0;
j<
j++)
a[i][j];
请输入方程组右端的常数项:
cin>
b[i];
for(j=0;
N-1;
{p=j;
for(i=j+1;
i++)
//寻找系数矩阵每列的最大值
{if(fabs(a[i][j])>
fabs(a[p][j]))
p=i;
}
if(p!
=j)
//交换第p行与第j行
{for(k=0;
k<
k++)
{
t=a[j][k];
a[j][k]=a[p][k];
a[p][k]=t;
}
//交换常数项的第p行与第j行
t=b[p];
b[p]=b[j];
b[j]=t;
//把系数矩阵第j列a[j][j]下面的元素变为0
{
m=-a[i][j]/a[j][j];
for(n=j;
n<
n++)
a[i][n]=a[i][n]+a[j][n]*m;
b[i]=b[i]+b[j]*m;
}
//求方程组的一个解
x[N-1]=b[N-1]/a[N-1][N-1];
//回代法求方程组其他解
for(i=N-2;
i>
=0;
i--)
s=0.0;
for(j=N-1;
j>
i;
j--)
s=s+a[i][j]*x[j];
x[i]=(b[i]-s)/a[i][i];
方程组的解如下:
for(i=0;
=N-1;
x[i]<
2.4高斯列主元法解线性方程组程序调试结果图示:
求解下列方程组
图2-2高斯列主元法解线性方程组程序算法调试图
3.牛顿法解非线性方程组
3.1牛顿法解非线性方程组算法说明
牛顿法主要思想是用
进行迭代。
因此首先需要算出
的雅可比矩阵
,再求过
求出它的逆
,当它达到某个精度时即停止迭代。
具体算法如下:
首先设
已知。
①:
计算函数
(3-1)
②:
计算雅可比矩阵
(3-2)
(3-3)
③:
求线性方程组
的解
。
④:
计算下一点
重复上述过程。
3.2牛顿法解非线性方程组流程图
图3-1牛顿法解非线性方程组流程图
3.3牛顿法解非线性方程组程序代码
#defineN2
//非线性方程组中方程个数、未知量个数
#defineEpsilon0.0001
//差向量1范数的上限
#defineMax
100
//最大迭代次数
constintN2=2*N;
voidff(floatxx[N],floatyy[N]);
//计算向量函数的因变量向量yy[N]
voidffjacobian(floatxx[N],floatyy[N][N]);
//计算雅克比矩阵yy[N][N]
voidinv_jacobian(floatyy[N][N],floatinv[N][N]);
//计算雅克比矩阵的逆矩阵inv
voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N]);
//由近似解向量x0计算近似解向量x1
floatx0[N]={2.0,0.25},y0[N],jacobian[N][N],invjacobian[N][N],x1[N],errorno
rm;
inti,j,iter=0;
//如果取消对x0的初始化,撤销下面两行的注释符,就可以由键盘向x0读入初始近似解向量
//for(i=0;
//
x0[i];
初始近似解向量:
for(i=0;
x0[i]<
"
do
iter=iter+1;
第"
iter<
次迭代开始"
//计算向量函数的因变量向量y0
ff(x0,y0);
//计算雅克比矩阵jacobian
ffjacobian(x0,jacobian);
//计算雅克比矩阵的逆矩阵invjacobian
inv_jacobian(jacobian,invjacobian);
newdundiedai(x0,invjacobian,y0,x1);
//计算差向量的1范数errornorm
errornorm=0;
errornorm=errornorm+fabs(x1[i]-x0[i]);
if(errornorm<
Epsilon)break;
x0[i]=x1[i];
}while(iter<
Max);
return0;
voidff(floatxx[N],floatyy[N])
{floatx,y;
x=xx[0];
y=xx[1];
yy[0]=x*x-2*x-y+0.5;
yy[1]=x*x+4*y*y-4;
向量函数的因变量向量是:
for(i=0;
yy[i]<
voidffjacobian(floatxx[N],floatyy[N][N])
floatx,y;
inti,j;
//jacobianhaven*nelement
yy[0][0]=2*x-2;
yy[0][1]=-1;
yy[1][0]=2*x;
yy[1][1]=8*y;
雅克比矩阵是:
{for(j=0;
yy[i][j]<
voidinv_jacobian(floatyy[N][N],floatinv[N][N])
{floataug[N][N2],L;
inti,j,k;
开始计算雅克比矩阵的逆矩阵:
for(i=0;
for(j=0;
aug[i][j]=yy[i][j];
for(j=N;
N2;
if(j==i+N)aug[i][j]=1;
else
aug[i][j]=0;
aug[i][j]<
for(k=i+1;
{L=-aug[k][i]/aug[i][i];
for(j=i;
aug[k][j]=aug[k][j]+L*aug[i][j];
{
for(i=N-1;
0;
i--)
for(k=i-1;
k>
k--)
for(j=N2-1;
aug[i][j]=aug[i][j]/aug[i][i];
inv[i][j-N]=aug[i][j];
雅克比矩阵的逆矩阵:
inv[i][j]<
voidnewdundiedai(floatx0[N],floatinv[N][N],floaty0[N],floatx1[N])
floatsum=0;
{sum=0;
sum=sum+inv[i][j]*y0[j];
x1[i]=x0[i]-sum;
近似解向量:
x1[i]<
3.4牛顿法解非线性方程组程序调试结果图示:
图3-2牛顿法解非线性方程组程序算法调试图
图3-3牛顿法解非线性方程组程序算法调试图
图3-4牛顿法解非线性方程组程序算法调试图
4.龙贝格求积分算法
4.1龙贝格求积分算法分析
生成j>
=k的近似积分结果逼近表,并以R(j+1,j+1)为最终解来逼近积分。
R(j,0)=T(j),j>
=0,T(j)为区间逐次减半递推梯形求积分公式算出的结果;
R(j,1)=S(j),j>
=1,S(j)为区间逐次减半递推辛普森求积分公式算出的结果;
(4-1)
R(j,2)=B(j),j>
=2,B(j)为递推布尔求积分公式算出的结果;
(4-2)
生成
的逼近表
,并以
为最终解来逼近积分
(4-3)
逼近
存在于一个特别的下三角矩阵中,第0列元素
用基于
个[a,b]子区间的连续梯形方法计算,然后利用龙贝格公式计算
当
时,第
行的元素为
时,程序在第
行结束。
4.2龙贝格积分算法流程图
图4-1龙贝格积分算法流程图
4.3龙贝格积分算法程序代码
stdio.h>
math.h>
#definef(x)(sin(x))
//列举函数
#defineN_H20
#defineMAX
10
#define
a
0
//所求积分的上下限
b
1
epsilon
0.0001
//所需精度
doubleRomberg(doubleaa,doublebb,longintn)
inti;
doublesum,h=(bb-aa)/n;
sum=0;
for(i=1;
n;
sum+=f(aa+i*h);
sum+=(f(aa)+f(bb))/2;
return(h*sum);
longintn=N_H,m=0;
doubleT[2][MAX+1];
T[1][0]=Romberg(a,b,n);
n*=2;
for(m=1;
m<
MAX;
m++)
m;
{T[0][i]=T[1][i];
n=n*2;
for(i=1;
=m;
T[1][i]=T[1][i-1]+(T[1][i-1]-T[0][i-1])/(pow(2,2*m)-1);
if((T[0][m-1]-T[1][m])<
epsilon)
T="
T[1][m]<
return;
4.4龙贝格积分算法调试图
求
在
区间上的精度为0.0001的积分
图4-2龙贝格积分程序算法调试图
5.三次样条插值算法
5.1三次样条插值基本算法说明:
表5-1三次样条插值基本算法说明
策略描述
包含
和
的方程
(i)三次紧压样条,确定
,
(如果导数已知,这是“最佳选择”)
(ii)natural三次样条(一条“松弛曲线”)
(iii)外挂
到端点
若已知N+1个点的
及其一阶导数的边界条件S’(a)=
和S’(b)=
则存在唯一的三次样条曲线。
构造并求解下列线性方程组
(5-1)
(5-2)
当得到系数{
}后,可利用如下公式计算样条函数
(5-3)
为了更有效的计算,每个三次多项式可表示成嵌套乘的形式,也可以写为如下形式:
(5-4)
5.2三次样条插值算法程序代码
stdlib.h>
#defineMAX4
double*diff(doubleX[],intn)
inti=0;
d