完整word版西安交通大学计算方法A上机大作业word文档良心出品.docx
- 文档编号:812157
- 上传时间:2022-10-13
- 格式:DOCX
- 页数:20
- 大小:397.59KB
完整word版西安交通大学计算方法A上机大作业word文档良心出品.docx
《完整word版西安交通大学计算方法A上机大作业word文档良心出品.docx》由会员分享,可在线阅读,更多相关《完整word版西安交通大学计算方法A上机大作业word文档良心出品.docx(20页珍藏版)》请在冰豆网上搜索。
完整word版西安交通大学计算方法A上机大作业word文档良心出品
计算方法A上机大作业
张晓璐硕4011班学号:
3114009097
1.共轭梯度法求解线性方程组
算法原理:
由定理3.4.1可知系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求解二次函数
极小点具有等价性,所以可以利用共轭梯度法求解
的极小点来达到求解Ax=b的目的。
共轭梯度法在形式上具有迭代法的特征,在给定初始值情况下,根据迭代公式:
产生的迭代序列
在无舍入误差假定下,最多经过n次迭代,就可求得
的最小值,也就是方程Ax=b的解。
首先导出最佳步长
的计算式。
假设迭代点
和搜索方向
已经给定,便可以通过
的极小化
来求得,根据多元复合函数的求导法则得:
令
,得到:
,其中
然后确定搜索方向
。
给定初始向量
后,由于负梯度方向是函数下降最快的方向,故第一次迭代取搜索方向
。
令
其中
。
第二次迭代时,从
出发的搜索方向不再取
,而是选取
,使得
与
是关于矩阵A的共轭向量,由此可求得参数
:
然后从
出发,沿
进行搜索得到
设已经求出
,计算
。
令
,选取
,使得
和
是关于A的共轭向量,可得:
具体编程计算过程如下:
(1)给定初始近似向量
以及精度
;
(2)计算
,取
;
(3)Fork=0ton-1do
(i)
;
(ii)
;
(iii)
;
(iv)若
,则输出近似解
,停止;否则,转(v);
(v)
;
(vi)
;
Enddo
程序框图:
程序使用说明:
本共轭梯度法求解线性方程的程序是一个函数Gongetidu2(A,b,x0,epsa),在求解线性方程组Ax=b(A是对称正定矩阵)的时候,首先给定初始向量x0和误差epsa,然后直接调用该函数Gongetidu2(A,b,x0,epsa)即可,其中A,b,x0和epsa都是可变的,虽然该函数是通用的,但是对于矩阵A和向量b的输入,需要使用者根据A和b的特点自行输入。
算例3.2计算结果:
对题113页3.2,首先编程将矩阵A,b,x0,epsa读入系统,然后再调用函数Gongetidu2(A,b,x0,epsa);
程序如下:
clearAbx0%程序运行前首先清除A,b和X0的数值,以免影响计算
clc
n=input('请输入对称正定矩阵A的阶数n=');
epsa=input('请输入误差=');
fork=1:
(n-1)%读取题目3.2的矩阵A
A(k,k)=-2;
A(k,k+1)=1;
A(k+1,k)=1;
end
A(n,n)=-2;
A;
b(1,1)=-1;%读取题目3.2的矩阵b
b(n,1)=-1;
b;
x0(1,1)=input('假定初始向量各元素相同,均为:
');%给题目3.2的迭代过程赋初值
forkk=2:
n
x0(kk,1)=x0(1,1);
end
x=Gongetidu2(A,b,x0,epsa)%调用共轭梯度法求线性方程组的函数
functionx=Gongetidu2(A,b,x0,epsa)
n=size(A,1);
x=x0;
r=b-A*x;
d=r;
fork=0:
(n-1)
alpha=(r'*r)/(d'*A*d);
x=x+alpha*d;
r2=b-A*x;
if((norm(r2)<=epsa)|(k==n-1))
x;
break;
end
beta=norm(r2)^2/norm(r)^2;
d=r2+beta*d;
r=r2;
end
计算结果:
当n=100时,x=[1;1;1;1;1;1;1;…..1;1;1;1]
当n=200时,x=[1;1;1;1;1;1;1;…..1;1;1;1;1;1;1]
当n=400时,x=[1;1;1;1;1;1;1;…..1;1;1;1;1;1;1;1;1;1]
2.三次样条差值
算法原理(三次样条插值函数的导出):
(i).导出在子区间
上的S(x)的表达式
由于S(x)的二阶导数连续,设S(x)再节点
处的二阶导数值S’’(xi)=Mi,其中Mi为未知的待定参数。
由S(x)是分段的三次多项式知,S’’(x)是分段线性函数,S’’(x)在子区间
上可表示为
其中hi=xi-x(i-1),对上式两次积分得到
由插值条件
得到
将
代入
可得
(ii).建立参数
的方程组
对
S(x)求导可得
上式中令
得S(x)在xi处的左导数
,令
得到右导数
,因为S(x)在内节点xi处一阶导数连续,所以
,进一步推导可得
其中
,
,
上式为三弯矩方程组,因为三弯矩方程组只有n-1个方程,不能确定n+1个未知量Mi,所以需要再增加两个方程,由边界条件确定。
第一种边界条件:
此时已知
.不妨取
,
,这时三弯矩方程组化为:
以上方程组系数矩阵式严格三对角占优矩阵,可用追赶法求解。
求出
后,代入S(x)可得三次样条插值函数的数学表达式。
第二种边界条件:
已知
。
记
,则有
所以:
即
其中
所以得到第二种边界条件下的三弯矩方程组:
该方程组系数矩阵是严格三对角占优矩阵,可用追赶法求解,具体追赶法的求解过程见《数值分析》教材。
第三种边界条件:
周期型边界条件.已知
是以
为周期的周期函数,则由周期性可知,
,这时,可以将点
看成内点,则方程组对i=n也成立,既有
,
也即
,
其中
于是三弯矩方程组化为
该方程组可用LU分解法求解,具体追赶法的求解过程见《数值分析》教材。
程序框图如下:
程序使用说明:
因为该程序需要计算三种不同边界条件的三弯矩方程组,所以首先定义追赶法和LU分解法求解线性方程组的函数,ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b),然后在确定了边界条件类型之后,构造关于M的三弯矩方程AM=b的A和b矩阵,然后分别代入函数ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b)便可求得M,然后根据
便可得到,在
上的三次样条插值函数
,进而得到整个区间上的三次样条差值函数
。
算例计算结果:
141页的计算实习
当为第一类边界条件时:
当-1= 当-0.8= 当-0.6= 当-0.4= 当-0.2= 当0= 当0.2= 当0.4= 当0.6= 当0.8= 当为第二类边界条件时: 当-1= 当-0.8= 当-0.6= 当-0.4= 当-0.2= 当0= 当0.2= 当0.4= 当0.6= 当0.8= 3.龙贝格积分法: 对于复化梯形求积公式,取积分近似值 对复化辛普森求积公式 因为 ,所以上述两式相除得 所以, 同理, 对 , 和 分析可得 龙贝格积分算法如下: 如果 则取 ,否则,继续计算直到满足条件。 程序框图: 程序使用说明: 运行本程序的时候,直接按照提示输入所求积分的原函数 (比如 ),然后按提示依次输入积分下限 ,积分上限 和积分精度 ,然后程序便可计算出原函数 在 之间的积分数值。 算例计算结果: 6.2 4.四阶龙格-库塔法求解常微分方程的初值问题 算法原理: 用标准4级4阶R-K法求解一阶常微分方程的算法公式为: 求解m阶常微分方程 将其转化为一阶微分方程组来求解。 为此,引进新的变量 …, 即可将m阶微分方程,转化为如下的一阶微分方程组: 程序框图: 程序使用说明: 运行该程序,按照提示依次输入常微分方程的阶数m,x的下限a和上限b,步长h以及y的最高阶导数表达式 ,特别需要注意的是: 由于程序设定的问题,一定要将公式中的 用 代替,例如所求解的初值问题为 ,那么本程序所需要输入的公式为: 。 且 代表含 的表达式。 接着依次输入已知条件 , , ,…,即可,最终可以得出该问题的解。 算例计算结果: 9.2题结果
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 完整 word 西安交通大学 计算方法 上机 作业 文档 良心 出品
![提示](https://static.bdocx.com/images/bang_tan.gif)