计算方法
试验报告
班级:
学号:
姓名:
实验牛顿山法
1 实验目
(1) 熟悉非线性方程求根简单迭代法牛顿迭代牛顿山法
(2) 编程实现简单迭代法牛顿迭代牛顿山法
(3) 认识选择迭代格式重性
(4) 迭代速度建立感性认识分析实验结果体会初值迭代影响
2 实验容
(1)牛顿山法解方程(初值06)
输入:初值误差限迭代次数山次数
输出:似根步山子
(2)设方程f(x)x 3x –10 三实根 x18793 x034727 x153209现采面六种计算格式求 f(x)0根 x x
x x x
x x x x
输入:初值误差限迭代次数
输出:似根实际迭代次数
3 算法基原理
求非线性方程组解科学计算常遇问题实际背景.种算法层出穷中迭代流算法建立效迭代格式迭代数列收敛求根设计算法前般迭代进行收敛性判断关重牛顿法切线法迭代算法中典型方法初值选取适单根附牛顿法收敛速度快初值牛顿迭代
关重初值选取采牛顿山算法进行纠正
般迭代:
牛顿公式:
牛顿山公式:
图31般迭代算法流程图
山子
山条件
4 算法描述
般迭代算法见流程图
图32牛顿山算法流程图
牛顿山算法见流程图:
5代码:
#include
#include
#include
using namespace std
class srrt
{
private
int n
double *a *xr *xi
public
srrt (int nn)
{
n nn
a new double[n+1] 动态分配存
xr new double[n]
xi new double[n]
}
void input () 文件读入代数方程系数
void srrt_root () 执行牛顿山法
void output () 输出根文件显示
~srrt ()
{ delete [] a xr xi }
}
void srrtinput () 文件读入代数方程系数
{
int i
char str1[20]
cout <<\n输入文件名
cin >>str1
ifstream fin (str1)
if (fin)
{ cout <<\n开文件 <
finclose ()
}
void srrtsrrt_root () 执行牛顿山法
{
int mijtkisit
double txyx1y1dxdypqwdddcc
double guvpqg1u1v1
mn
while ((m>0)&&(fabs(a[m])+1010)) mm1
if (m<0)
{
cout <<\n程序工作失败 <
}
for (i0 i
wa[i] a[i]a[mi] a[mi]w
}
km is0 w10
jt1
while (jt1)
{
pqfabs(a[k])
while (pq<10e12)
{
xr[k1]00 xi[k1]00 kk1
if (k1)
{
xr[0]a[1]*wa[0] xi[0]00
return
}
pqfabs(a[k])
}
qlog(pq) qq(10*k) qexp(q)
pq ww*p
for (i1 i
x00001 x1x y02 y1y dx10
g10e+37
l40
ua[0] v00
for (i1 i
pu*x1 qv*y1
pq(u+v)*(x1+y1)
upq+a[i] vpqpq
}
g1u*u+v*v
if (g1>g)
{
if (is0)
{
it1
if (it0)
{
is1
ddsqrt(dx*dx+dy*dy)
if (dd>10) dd10
dc628(45*k) c00
}
while(11)
{
cc+dc
dxdd*cos(c) dydd*sin(c)
x1x+dx y1y+dy
if (c<629) { it0 break }
dddd167
if (dd<10e07) { it1 break }
c00
}
if (it0) goto l40
}
else
{
it1
while (it1)
{
tt167 it0
x1xt*dx
y1yt*dy
if (k>50)
{
psqrt(x1*x1+y1*y1)
qexp(850k)
if (p>q) it1
}
}
if (t>10e03) goto l40
if (g>10e18)
{
it0
if (it0)
{
is1
ddsqrt(dx*dx+dy*dy)
if (dd>10) dd10
dc628(45*k) c00
}
while(11)
{
cc+dc
dxdd*cos(c) dydd*sin(c)
x1x+dx y1y+dy
if (c<629) { it0 break }
dddd167
if (dd<10e07) { it1 break }
c00
}
if (it0) goto l40
}
}
if (fabs(y)<10e06)
{ px y00 q00 }
else
{
p20*x qx*x+y*y
xr[k1]x*w
xi[k1]y*w
kk1
}
for (i1 i
a[i]a[i]a[i1]*p
a[i+1]a[i+1]a[i1]*q
}
xr[k1]x*w xi[k1]y*w
kk1
if (k1)
{ xr[0]a[1]*wa[0] xi[0]00 }
}
else
{
gg1 xx1 yy1 is0
if (g<10e22)
{
if (fabs(y)<10e06)
{ px y00 q00 }
else
{
p20*x qx*x+y*y
xr[k1]x*w
xi[k1]y*w
kk1
}
for (i1 i
a[i]a[i]a[i1]*p
a[i+1]a[i+1]a[i1]*q
}
xr[k1]x*w xi[k1]y*w
kk1
if (k1)
{ xr[0]a[1]*wa[0] xi[0]00 }
}
else
{
u1k*a[0] v100
for (i2 i
pu1*x qv1*y pq(u1+v1)*(x+y)
u1pq+(ki+1)*a[i1]
v1pqpq
}
pu1*u1+v1*v1
if (p<10e20)
{
it0
if (it0)
{
is1
ddsqrt(dx*dx+dy*dy)
if (dd>10) dd10
dc628(45*k) c00
}
while(11)
{
cc+dc
dxdd*cos(c) dydd*sin(c)
x1x+dx y1y+dy
if (c<629) { it0 break }
dddd167
if (dd<10e07) { it1 break }
c00
}
if (it0) goto l40
if (fabs(y)<10e06)
{ px y00 q00 }
else
{
p20*x qx*x+y*y
xr[k1]x*w
xi[k1]y*w
kk1
}
for (i1 i
a[i]a[i]a[i1]*p
a[i+1]a[i+1]a[i1]*q
}
xr[k1]x*w xi[k1]y*w
kk1
if (k1)
{ xr[0]a[1]*wa[0] xi[0]00 }
}
else
{
dx(u*u1+v*v1)p
dy(u1*vv1*u)p
t10+40k
it1
while (it1)
{
tt167 it0
x1xt*dx
y1yt*dy
if (k>50)
{
psqrt(x1*x1+y1*y1)
qexp(850k)
if (p>q) it1
}
}
if (t>10e03) goto l40
if (g>10e18)
{
it0
if (it0)
{
is1
ddsqrt(dx*dx+dy*dy)
if (dd>10) dd10
dc628(45*k) c00
}
while(11)
{
cc+dc
dxdd*cos(c) dydd*sin(c)
x1x+dx y1y+dy
if (c<629) { it0 break }
dddd167
if (dd<10e07) { it1 break }
c00
}
if (it0) goto l40
}
if (fabs(y)<10e06)
{ px y00 q00 }
else
{
p20*x qx*x+y*y
xr[k1]x*w
xi[k1]y*w
kk1
}
for (i1 i
a[i]a[i]a[i1]*p
a[i+1]a[i+1]a[i1]*q
}
xr[k1]x*w xi[k1]y*w
kk1
if (k1)
{ xr[0]a[1]*wa[0] xi[0]00 }
}
}
}
if (k1) jt0
else jt1
}
}
void srrtoutput () 输出根文件显示
{
int k
char str2[20]
cout <<\n输出文件名
cin >>str2
ofstream fout (str2)
if (fout)
{ cout <<\n开文件 <
fout <
foutclose ()
}
void main () 函数
{
srrt root(6)
rootinput () 文件读入代数方程系数
rootsrrt_root () 执行牛顿山法
rootoutput () 输出根文件显示
}
6输入输出
输出结果:
7分析体会
牛顿山法作计算方法课程中重知识点书学时较易致理解思路级编写代码时难度学鉴代码够实现基功实验程中加深牛顿山法理解
实验二 高斯消法
1 实验目
(1) 熟悉求解线性方程组关理方法
(2) 编程实现雅高斯塞德尔迭代法列元高斯消法约消追赶法
(3) 通测试进步解种方法优缺点
(4) 根类型方程组选择合适数值方法
2 实验容
(1)Gauss Seidel 迭代法求解方程组
输入:系数矩阵A迭代次数N初始量误差限e
输出:解量
(2)城镇三生产企业:煤矿电厂方铁路作济系统
煤矿:生产价值1元煤需消耗025元电费035元运输费
电厂:生产价值1元电需消耗040元煤费005元电费010元运输费
铁路:提供价值1元铁路运输服务需消耗045元煤010元电费010元运输费
某星期三企业间彼需求煤矿50000元订单电厂25000元电量供应求方铁路价值30000元运输需求
提问:三企业星期应生产少产值满足外需求?
提示:
投入产出间衡关系:物资消耗新创造价值等生产总产值
输入:物资消耗新创造价值
输出:三企业生产总值
(3)天文学家确定颗行星绕太阳运行轨道轨道面建立太阳原点直角坐标系两坐标轴取天文测量单位(1天文单位球太阳均距离:9300万里)五时间行星作五次观测轨道五点坐标分(57640648)(62861202)(67591823)(71682562)(74083360)开普勒第定律知行星轨道椭圆试建立方程
(4)选元高斯消求行列式值
提示:
A
B 消元结果直接存储系数矩阵中
C 消元程发生两行调情况偶数次时行列式值角线积否角线积相反数
(5)选元约消分矩阵求逆矩阵逆输出奇异标志
提示:
3 算法基原理
三次样条拟合问题终结线性方程组求解线性方程组数值分析中非常重工程计算中容忽视
线性方程组致分迭代法直接法收敛条件满足时进行迭代雅高斯塞德尔基两类迭代方法区迭代程中否引新值进行剩计算消元简单直接法十分效列元高斯消法求解般线性方程组适时求矩阵应行列式约消实质初等行变换系数矩阵化单位阵求矩阵逆直接法注意空间时间两方面算法进行优化
高斯塞德尔迭代:
列元高斯消法:列元
消元 回代
图42列元约消
图41GS迭代算法流程图
约消
4 算法描述
Gauss Seidel算法见流程图
选列元高斯消法见流程图
选列元约消法见流程图
5计算例参考输出
实验1输出参考书p146页 表52
实验2参考输出:三企业产值X[114458 65395 85111]
实验3参考输出:椭圆曲线系数a[00508 00702 00381 04531 02643]
实验4参考输出:18
实验5参考输出:
图43列元高斯消流程图
逆
6代码
#include
#include
#define N 3
using namespace std
int main()
{
ifstream file(datxt)
double x[N][N+1]
double result[N]
for(int i0i
for(int j0j
file>>x[i][j]
}
}
for(i0i
for(int ji+1j
x[i][j]x[i][j]x[i][i]
}
for(ji+1j
for(int mi+1m
x[j][m]x[j][m]x[j][i]*x[i][m]
}
}
}
for(iN1i>0i)
{
result[i]x[i][N]
for(int jN1j>ij)
{
result[i]result[i]result[j]*x[i][j]
}
}
for(i0i
}
7输出
实验结果:
8分析体会
消元程中次选择绝值者作元高斯消法重点高斯消法计算高次较快捷方便种算法解决线性方程组方程组化三角形格式高斯消元法找出解学高斯消法解线性方程组十分重
实验三 高斯赛德尔迭代
1 实验目
(5) 熟悉求解线性方程组关理方法
(6) 编程实现雅高斯塞德尔迭代法列元高斯消法约消追赶法
(7) 通测试进步解种方法优缺点
(8) 根类型方程组选择合适数值方法
2 实验容
(1)Gauss Seidel 迭代法求解方程组
输入:系数矩阵A迭代次数N初始量误差限e
输出:解量
(2)城镇三生产企业:煤矿电厂方铁路作济系统
煤矿:生产价值1元煤需消耗025元电费035元运输费
电厂:生产价值1元电需消耗040元煤费005元电费010元运输费
铁路:提供价值1元铁路运输服务需消耗045元煤010元电费010元运输费
某星期三企业间彼需求煤矿50000元订单电厂25000元电量供应求方铁路价值30000元运输需求
提问:三企业星期应生产少产值满足外需求?
提示:
投入产出间衡关系:物资消耗新创造价值等生产总产值
输入:物资消耗新创造价值
输出:三企业生产总值
(3)天文学家确定颗行星绕太阳运行轨道轨道面建立太阳原点直角坐标系两坐标轴取天文测量单位(1天文单位球太阳均距离:9300万里)五时间行星作五次观测轨道五点坐标分(57640648)(62861202)(67591823)(71682562)(74083360)开普勒第定律知行星轨道椭圆试建立方程
(4)选元高斯消求行列式值
提示:
D
E 消元结果直接存储系数矩阵中
F 消元程发生两行调情况偶数次时行列式值角线积否角线积相反数
(5)选元约消分矩阵求逆矩阵逆输出奇异标志
提示:
3 算法基原理
三次样条拟合问题终结线性方程组求解线性方程组数值分析中非常重工程计算中容忽视
线性方程组致分迭代法直接法收敛条件满足时进行迭代雅高斯塞德尔基两类迭代方法区迭代程中否引新值进行剩计算消元简单直接法十分效列元高斯消法求解般线性方程组适时求矩阵应行列式约消实质初等行变换系数矩阵化单位阵求矩阵逆直接法注意空间时间两方面算法进行优化
高斯塞德尔迭代:
列元高斯消法:列元
消元 回代
图42列元约消
图41GS迭代算法流程图
约消
4 算法描述
Gauss Seidel算法见流程图
选列元高斯消法见流程图
选列元约消法见流程图
5计算例参考输出
实验1输出参考书p146页 表52
实验2参考输出:三企业产值X[114458 65395 85111]
实验3参考输出:椭圆曲线系数a[00508 00702 00381 04531 02643]
实验4参考输出:18
实验5参考输出:
图43列元高斯消流程图
逆
6代码:
#include
#include
using namespace std
int main(){
double xishu[3][4]
double x[3]{0}
ifstream in(D1txt)
for(int i0i<3i++)
{
for(int j0j<4j++)
{
in>>xishu[i][j]
}
}
inclose()
ofstream out(d2txt)
int t0
while(t++<7){
for(int i0i<3i++){
x[i]xishu[i][3]
for(int j0j<3j++){
if((ji))
{
x[i]xishu[i][j]*x[j]
}
}
x[i]xishu[i][i]
out<
out<
return 0
}
7输入输出
8分析体会
次试验通理实例线性方程组解法收敛性误差分析进行探讨线性方程组数值解法讨高斯赛德尔迭代法进步研究总结高斯赛德尔迭代法理应分析问题编辑程序时更握高斯赛德尔迭代法应
实验四 二法三次样条插值
1 实验目
(1) 明确插值项式分段插值项式优缺点
(2) 编程实现拉格朗日插值算法分析实验结果体会高次插值产生龙格现象
(3) 理解二拟合编程实现线性拟合掌握非线性拟合转化线性拟合方法
(4) 运常插值拟合方法解决实际问题
2 实验容
(1)
求选取11等距插值节点分采拉格朗日插值分段线性插值计算x05 45处函数值结果精确值进行较
输入:区间长度n(n+1节点)预测点
输出:预测点似函数值精确值误差
(2)定元函数 n+1节点值
数:
04
055
065
080
095
105
041075
057815
069675
090
100
125382
求五次拉格朗日项式分段三次插值项式计算 值
输入:数点集数预测点
输出:预测点似函数值
(3)某冶炼程中根统计数含碳量时间关系试求含碳量时间拟合曲线
分)
0 5 10 15 20 25 30 35 40 45 50 55
0 127 216 286 344 387 415 437 451 458 402 464
输入: 数点集数
输出:拟合函数印出误差
提示:似解析表达式
(4)线性拟合编程实现药方案
背景知识:种新药床前必须设计药方案
床种药物效浓度c1效浓度c2
设计药方案时血药浓度 保持c1~c2间 题设c110c225(ugml)
实验方面某快速静脉注射方式次注入该药物300mg定时刻t(时)采集血药测血药浓度c(ugml)表
t (h) 025 05 1 15 2 3 4 6 8
c (mgml) 1921 1815 1536 1410 1289 932 745 524 301
问题:
B 快速静脉注射药方式研究血药浓度(单位体积血液中药物含量)变化规律
C 设计药方案:次注射剂量间隔时间长
提示:
A 血药浓度变化规律符合指数关系
B 指数关系转化线性关系
3 算法基原理
精确函数 y f(x) 非常复杂未知系列节点 x0 … xn 处测函数值 y0 f(x0) … yn f(xn) 希构造简单易算似函数 g(x) » f(x)满足条件g(xi) f(xi) (i 0 … n)里 g(x) 称f(x) 插值函数插值函数似估计f(x)未知点处函数值常插值函数项式插值
谓项式插值求 n 次项式
基基函数拉格朗日插值构造项式插值基方法推导数值微积分微分方程数值解公式理基础
拉格朗日插值公式: 中
结点较次数较高插值项式发生Runge现象分段低次插值避免Runge现象重手段分段次插值整区间分段区间次项式逼 f (x)直观折线代曲线区间足够保证逼效果曲线缺乏光滑性
插值逼复杂函数f(x)种方法拟合结点数较者够精确采拟合方法拟合简单函数p(x)需严格通定结点结点处总体误差达极满足二拟合条件线性拟合简单拟合方法许非线性问题转化线性拟合令线性函数a0+a1x确定a0a1终求解线性方程组:
图11拉格朗日插值算法
4 算法描述
(2) 拉格朗日插值算法见流程图
(3) 分段线性插值算法注意关键点:
A 定位:预测点x第区间?
B 假定x第k区间
次插值公式
(3) 线性拟合算法
没讨线性方程组数值解法前提原始消元法求解线性方程组
5 计算例参考输出
(1)实验1参考输出
X y(精确) y(拉格朗日) y(分段线性) 误差(拉) 误差(分)
0500000 0800000 0843408 0750000 0043408 0050000
4500000 0047059 1578721 0048643 1537662 0001584
(2)实验2参考输出
模型中参数 k02347 v1502
制定药方案:
首次注射 375 mg 余次注射 225 mg 注射间隔时间 4 时
6 代码
二法:
#include
#include
using namespace std
void main()
{
ifstream in(d2txt)
float **axresult
int n
float kb
void smallest(float **aint nfloat*kfloat *b)
cout<<输入点数:
in>>n
anew float *[n]
for(int i0i
cout<<输入点:<
in>>a[i][j]
smallest(an&b&k)
cout<<输入变量值:
in>>x
resultb+k*x
out<<拟合曲线方程\nY< out<
void smallest(float **aint nfloat*kfloat *b)
{
float can[2][3]{n00000}
for(int i0i
can[0][1]+a[i][0]
can[0][2]+a[i][1]
can[1][1]+a[i][0]*a[i][0]
can[1][2]+a[i][0]*a[i][1]
}
can[1][0]can[0][1]
for(i2i>0i)
can[1][i]can[1][i]can[0][i]*can[1][0]can[0][0]
*bcan[1][2]can[1][1]
*k(can[0][2]can[0][1]**b)can[0][0]
}
三次样条插值:
#include
#include
#include
using namespace std
const int MAX 50
float x[MAX] y[MAX] h[MAX]
float c[MAX] a[MAX] fxym[MAX]
float f(int x1 int x2 int x3){
float a (y[x3] y[x2]) (x[x3] x[x2])
float b (y[x2] y[x1]) (x[x2] x[x1])
return (a b)(x[x3] x[x1])
} 求差分
void cal_m(int n){ 追赶法求解出弯矩量M……
float B[MAX]
B[0] c[0] 2
for(int i 1 i < n i++)
B[i] c[i] (2 a[i]*B[i1])
fxym[0] fxym[0] 2
for(i 1 i < n i++)
fxym[i] (fxym[i] a[i]*fxym[i1]) (2 a[i]*B[i1])
for(i n1 i > 0 i)
fxym[i] fxym[i] B[i]*fxym[i+1]
}
void printout(int n)
int main(){
int times0
int ni char ch
ifstream in(d\\1txt)
do{
cout<
for(i 0 i < n i++){
cout<
for(i 0 i < n i++) 求 步长
h[i] x[i+1] x[i]
cout<
float f0 f1
in>>t
switch(t){
case 1cout<
c[0] 1 a[n] 1
fxym[0] 6*((y[1] y[0]) (x[1] x[0]) f0) h[0]
fxym[n] 6*(f1 (y[n] y[n1]) (x[n] x[n1])) h[n1]
break
case 2cout<
c[0] a[n] 0
fxym[0] 2*f0 fxym[n] 2*f1
break
defaultcout<<\n定
}switch
for(i 1 i < n i++)
fxym[i] 6 * f(i1 i i+1)
for(i 1 i < n i++){
a[i] h[i1] (h[i] + h[i1])
c[i] 1 a[i]
}
a[n] h[n1] (h[n1] + h[n])
cal_m(n)
if(times0)
cout<<已知两端阶导数时:\n
else
cout<<已知两端二阶导数时:\n
printout(n)
times++
cout<
}while(ch 'y' || ch 'Y')
return 0
}
void printout(int n){
cout<
cout< *
cout<
if(t > 0)cout<
if(t > 0)cout<< + <
t (y[i] fxym[i]*h[i]*h[i]6)h[i]
if(t > 0)cout<<+ <
if(t > 0)cout<< + <
cout<
7实验运行结果:
二法:
三次样条差值:
8分析体会
实际应中时需解决问题中运算复杂运算量时需助计算机帮助解决实际问题利二法样条插值C++语言中编程实现帮助更理解计算方法级实践益进步学现实生活中需通已数发掘事物身规律者模拟出相应数学模型解决需述相关知识完成
实验五龙贝格算法
1 实验目
(1) 熟悉复化梯形方法复化Simpson方法梯形递推算法龙贝格算法
(2) 编程实现复化梯形方法复化Simpson方法梯形递推算法龙贝格算法
(3) 理解掌握适应算法收敛加速算法基思想
(4) 分析实验结果体会种方法精确度建立计算机求解定积分问题感性认识
2 实验容
(1) 龙贝格算法计算
输入:积分区间误差限
输出:序列TnSnCnRn积分结果(参考书P81表24)
(2)设计方案计算国土面积
计算瑞士国土面积首先图作测量西东方x轴南北方y轴选择方便原点西边界东边界x轴区间适划分干
段分点y方测出南边界点北边界点y坐标数表(单位mm)
x 70 105 130 175 340 405 445 480 560
y1 44 45 47 50 50 38 30 30 34
y2 44 59 70 72 93 100 110 110 110
x 610 685 765 805 910 960 1010 1040 1065
y1 36 34 41 45 46 43 37 33 28
y2 117 118 116 118 118 121 124 121 121
x 1115 1180 1235 1365 1420 1460 1500 1570 1580
y1 32 65 55 54 52 50 66 66 68
y2 121 122 116 83 81 82 86 85 68
瑞士图外形图(例尺18mm40km)
问题:
试测量数计算瑞士国土似面积精确值41288方公里较
提示:
A 计算国土面积问题转化求积分问题
B 积函数表格形式(仅提供组数点函数形式未知)梯形公式简单
3 算法基原理
许实际问题中常常需计算定积分值根微积分学基定理积函数f(x)区间[ab]连续找f(x)原函数F(x)便利牛顿莱布尼兹公式求积分值
实际中遇困难牛顿莱布尼兹公式
(1) 找初等函数表示原函数
(2) 然找原函数表达式复杂便计算
(3) f(x)测量计算表格函数
种种困难必研究积分数值计算问题
利插值项式 积分转化显然易算称插值型求积公式简单插值型求积公式梯形公式Simpson公式求积结点提供较分段少结点梯形公式Simpson公式称复化梯形公式复化Simpson公式步长未知通误差限控制区间逐次分半策略动选取步长方法称适应算法梯形递推公式出区间分半前递推关系梯形递推公式求梯形序列相邻序列值作线性组合Simpson序列 Simpson序列作线性组合柯特斯序列 柯特斯序列作线性组合龙贝格序列|R2R1|
图21龙贝格算法数加工流程
复化梯形公式
复化Simpson公式
梯形递推公式
加权均公式:
龙贝格算法加快误差收敛速度梯形序列O(h2) 提高龙贝格序列O(h8)
4 算法描述
(1) 梯形递推算法见流程图
图22梯形递推算法流程图
图23龙贝格算法流程图
(2) 龙贝格算法见流程图
5 计算例参考输出
实验1参考输出
图24 龙贝格参考输出
6代码
hba
T2T1(F(a)+F(b))*h2
进行梯形公式数装载
for(int i0i
T[i]T2
sum0
xa+h2
while(x {
sumF(x)+sum
x+h
}
T2(T1+h*sum)2
hh2
T1T2
}
进行梯形公式辛甫生转化
for(int i0i
进行辛甫生柯特斯转化
for(int i0i
进行柯特斯龙贝格转化
for(int i0i
7实验运行图:
8分析体会
龙贝格积分法进行似求积原理埃特金插值类似进行线性整合结果具高精度求积效果实际程中评判合理步长困难常采取变步长办法进行计算结果满足精度求龙贝格运算程:梯形公式装载数→辛甫生化→柯特斯化→龙贝格高精度化通实验牢记龙贝格算法演算程更体会方法高效高精度
实验六 欧拉算法龙格库塔亚姆斯
1 实验目
(1) 熟悉数值微分中Euler法改进Euler法RungKutta方法
(2) 编程实现Euler法改进Euler法RungKutta方法
(3) 通实验结果分析算法优缺点
(4) 明确步长算法影响理解变步长RungKutta方法
2 实验容
(1) 0 < x<1
取h01时Euler法改进Euler法RungKutta方法求数值解精确解进行较
输入:求解区间初值数值解数
输出:数值解
(2)型火箭初始质量900千克中包括600千克燃料火箭竖直发射时燃料15千克秒速率燃烧掉产生30000牛顿恒定推力.燃料时引擎关闭设火箭升整程中空气阻力速度方成正例系数04(千克米).重力加速度取98米秒2
输入:相关数
输出:引擎关闭瞬间火箭高度速度加速度火箭达高点时间高度.
提示:火箭飞行程方程:
初始条件
3 算法基原理
许科学技术问题中建立模型常常常微分方程形式表示然少数特殊类型常微分方程解析方法求精确解外出般方程解析解表达式困难似方法求数值解实际工作中常计算机求常微分方程数值解谓常微分方程数值解常微分方程初值问题计算出系列节点 a x0< x1<…< xn b 处未知函数 y(x)似值y0y1…yn找系列离散点(x0y0)(x1y1)(x2y2)…(xnyn)似满足常微分方程数值解法基思想差商代导数实现连续问题离散化选取差商代导数公式
改进欧拉公式常方法包括预测校正两步先欧拉公式进行预报预报值代入梯形公式进行校正达二阶精度通龙格库塔法获更高精度典龙格库塔法区间[xnxn+1]取四点四点斜率进行加权均作均斜率通泰勒公式寻找局部截断误差O(h5)(4阶精度)参数满足条件
改进欧拉公式: 预测
校正
四阶(典)龙格库塔公式
4 算法描述
(1) 改进欧拉算法见流程图
(2) 典龙格库塔算法见流程图
5 计算例参考输出
(1)实验1参考输出
图53典龙格库塔算法
图52改进欧拉算法
(2)实验2参考输出
关闭引擎时(秒)时火箭高度h83224米速度v2541728米秒加速度a 40616米秒火箭高点时间51秒高度91864米
6 代码
改进欧拉格式:
*取h02改进欧拉方法求解列初值问题
dydxsqrt(2*x*x+3*y*y)y(0)50
#include
#include
#define MAXSIZE 50
double f(double xdouble y)
void main (void)
{
double abhx[MAXSIZE]y[MAXSIZE]
long in
printf(\n请输入求解区间ab)
scanf(1f1f&a&b)
printf(\n请输入步长h)
scanf(1f&h)
n(long)((ba)h)
x[0]a
printf(\n 请输入起点x[0]1f处坐标y[0]x[0])
scanf(1f&y[0])
for(i0i
x[i+1]x[i]+h
y[i+1]y[i]+h*f(x[i]y[i])
y[i+1]y[i]+h*(f(x[i]y[i])+f(x[i+1]y[i+1]))2
}
printf(\n 计算结果:)
for(i0i
}
double f(double xdouble y)
{
return(sqrt(2*x*x+3*y*y))
}
龙哥库塔:
#include
using namespace std
main()
{
int k
float X[k]Y[k]F[k]
float h0001
double K1K2K3K4
float f(floatfloat)
X[0]00
Y[0]10
F[0]f(X[0]Y[0])
for(k0k<100k++)
{ X[k+1]X[k]+h
K1h*F[k]
K2h*f(X[k]+05*hY[k]+05*K1)
K3h*f(X[k]+05*hY[k]+K2)
K4h*f(X[k]+hY[k]+K3)
Y[k+1]Y[k]+(K1+20*K2+20*K3+K4)60
F[k+1]f(X[k+1]Y[k+1])
cout<
system(Pause)
return 0
}
float f(float afloat b)
{ float c
cba*20b
return (c)
}
7
欧拉算法
龙格库塔
亚姆斯
8分析体会
改进欧拉算法龙格库塔算法亚姆斯算法计算方法中难点机实践掌握起会更透彻点程序编制程中体会抽象数值计算方法变具体计算机程序易程种数值计算方法数学公式更进步认识编制出计算机程序必须清楚数值计算方法具体计算步骤
学门课感觉实性较门课程连接数学计算机间桥梁计算方法门学科计算机结合体现出价值助计算机快速解决种方程求值结果具高度精确性单数学学做
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容,不代表《香当网》观点或立场,请自行判断内容的真实性和可靠性!
该内容是文档的文本内容,更好的格式请下载文档