数值分析 课程设计Word文档格式.docx
- 文档编号:19491173
- 上传时间:2023-01-06
- 格式:DOCX
- 页数:21
- 大小:196.14KB
数值分析 课程设计Word文档格式.docx
《数值分析 课程设计Word文档格式.docx》由会员分享,可在线阅读,更多相关《数值分析 课程设计Word文档格式.docx(21页珍藏版)》请在冰豆网上搜索。
。
(3)如果
,则停止计算(此时det(A)=0);
否则进行回代,
对于k=n,n-1,…,1,置
(当下标大于上标时,不做求和运算)。
(4)置
(5)输出
为线性方程组Ax=b的解,det为矩阵A的行列式值。
function[x,det,index]=Gauss(A,b)
%求线性方程组的列主元Gauss消去法,其中,
%A为方程组的系数矩阵;
%b为方程组的右端项;
%x为方程组的解;
%det为系数矩阵A的行列式的值;
%index为指标变量,index=0表示计算失败,index=1表示计算成功。
[n,m]=size(A);
nb=length(b);
%当方程组行与列的维数不相等时,停止计算,并输出出错信息。
ifn~=m
error('
TherowsandcolumnsofmatrixAmustbeequal!
'
);
return;
end
%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息。
ifm~=nb
ThecolumnsofAmustbeequalthelengthofb!
%开始计算,先赋初值
index=1;
det=1;
x=zeros(n,1);
fork=1:
n-1
%选主元
a_max=0;
fori=k:
n
ifabs(A(i,k))>
a_max
a_max=abs(A(i,k));
r=i;
end
ifa_max<
1e-10
index=0;
%交换两行
ifr>
k
forj=k:
z=A(k,j);
A(k,j)=A(r,j);
A(r,j)=z;
z=b(k);
b(k)=b(r);
b(r)=z;
%消元过程
fori=k+1:
m=A(i,k)/A(k,k);
forj=k+1:
A(i,j)=A(i,j)-m*A(k,j);
b(i)=b(i)-m*b(k);
det=det*A(k,k);
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<
index=0;
fork=n:
-1:
1
b(k)=b(k)-A(k,j)*x(j);
x(k)=b(k)/A(k,k);
实例1:
用列主元Gauss消去法解方程组
解:
在MATLAB命令窗口输入:
>
A=[-326;
10-70;
5-15];
b=[476]'
;
[x,det,index]=gauss(A,b)
x=
0.0000
-1.0000
1.0000
det=
155
index=
1
方程组的解为x=(0,-1,1)’,行列式值为det(A)=155,index=1表明计算是成功的。
实例2:
A=[10787;
7565;
86109;
75910];
b=[32233331]'
[x,det,index]=gauss(A,b)
运行结果:
方程组的解为x=(1,1,1,1)’,行列式值为det(A)=1,index=1表明计算是成功的。
主元素法有两种,一种列主元法,一种全主元法,一般来说列主元法就能确保算法的稳定,所谓算法的稳定是指在运算过程中计算误差(对消去法这种直接法来说主要指由于计算机字长有限带来的舍入误差)能得到控制,全主元是较列主元法更稳定的算法,但它的计算量要比列主元法大的多,列主元法在每做一次消元仅与同列的元素做比较,比较的次数与线性方程组的阶n是同阶的量,而全主元法每做一次消元要与系数矩阵所有元素进行比较,计算量是与n^2同阶的量,计算量较列主元大的多,一般来说不采用全主元法,而采用列主元法即可.
实验二Newton插值
1.熟悉Newton插值法的基本原理和方法
2.通过实例掌握用MATLAB求插值的方法
3.能编程实现Newton插值
根据均差定义,把x看成[a,b]上一点,可得
只要把后一式代入前一式,就得到
其中
是由上式定义的.确定的多项式
显然满足插值条件,且次数不超过n,它就是多项式
称为牛顿基本插值多项式,常用来计算非等距节点上的函数值。
MATLAB程序
functionN=newton(x,y,xi)
%Newton基本插值方法
n=length(x);
m=length(y);
ifm~=n
error('
xory输入有误,再来'
A=zeros(m);
Z=1.0;
A(:
1)=y;
N=A(1,1);
fork=2:
n%k为列标
fori=k:
n%i为行标
A(i,k)=(A(i-1,k-1)-A(i,k-1))/(x(i+1-k)-x(i));
Z=Z*(xi-x(k-1));
N=N+Z*A(k,k);
disp('
差商表:
一阶差商,二阶差商,三阶差商'
disp(A);
%disp(N);
fprintf('
Newton插值的结果保留3位小数是%10.3f\n'
N);
实验内容及结论
实例:
已知
,试用牛顿基本插值多项式计算
的近似值
x=[100121144]'
y=[101112]'
xi=151;
N=newton(x,y,xi)
差商表:
一阶差商,二阶差商,三阶差商
10.000000
11.00000.04760
12.00000.0435-0.0001
Newton插值的结果保留3位小数是12.285
N=12.2846
在Newton插值中运用牛顿基本插值多项式,并引入差商进行计算,可以用来计算非等距节点上的函数值。
实验三龙贝格求积公式
a.熟悉数值求积的龙本格求积法基本原理
b.掌握数值求积的龙本格求积算法程序设计
(1)置N=1,精度要求
,h1=b-a。
(2)用梯形公式计算积分近似值
(3)按变步长梯形公式计算积分近似值,将区间逐次分半,置
,并计算
(4)置M=N,N=2N,k=1。
(5)按加速公式求加速值:
(当k=1时被称为梯形加速公式,当k=2时,被称为辛普森加速公式,当k=3时,被称为龙贝格求积公式)。
(6)若M=1,转(7);
否则,置
,k=k+1,转(5)。
(7)若
,则停止计算(输出
),否则转(3)。
functionI=R_quad_iter(fun,a,b,ep)
%Romberg求积公式,其中,
%fun为被积函数;
%a,b为积分区间的端点,要求a<
b;
%ep为精度要求,默认值1e-5。
ifnargin<
4ep=1e-5;
m=1;
h=b-a;
I=h/2*(feval(fun,a)+feval(fun,b));
%用梯形公式计算积分近似值
T(1,1)=I;
while1
N=2^(m-1);
h=h/2;
I=I/2;
fori=1:
N
I=I+h*feval(fun,a+(2*i-1)*h)%按变步长梯形公式计算积分近似值将区间逐次分半,令区间长度为h/2,计算
T(m+1,1)=I;
M=2*N;
k=1;
whileM>
T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);
M=M/2;
k=k+1;
ifabs(T(k,k)-T(k-1,k-1))<
ep
break;
m=m+1;
I=T(k,k);
用Romberg求积法计算积分
,取精度要求
=10-5。
输入
fun=inline('
4/(1+x^2)'
%被积函数=内联函数
I=R_quad_iter(fun,0,1)
I=3.1000
I=2.4912
I=3.1312
I=2.0579
I=2.4963
I=2.8558
I=3.1390
I=1.8185
I=2.0600
I=2.2878
I=2.4976
I=2.6875
I=2.8573
I=3.0079
I=3.1409
I=3.1416
故最终结果为I=3.1416
在本积分公式中,将收敛缓慢的梯形值序列加工成收敛迅速的龙贝格值序列,使其精度更能满足实际要求。
实验四常微分方程数值解龙格—库塔
1.熟悉Runge-Kutta常微分方程初值问题的基本原理
2.了解Runge-Kutta常微分方程初值问题的计算流程
3.能编程实现Runge-Kutta常微分方程初值问题
常微分方程的解析解,由于求解的难度较高;
有时解析解太复杂,难以使用;
也有时的根无法取得解析解,所以在工程上往往求助于数值解。
所谓数值解,就是用加减乘除和函数求解等运算,在特定点求近似解的过程。
求解常微分方程的数值解的方法有欧拉法,改进的欧拉法和龙格—库塔法。
(1)欧拉法最为简单,它属于单步递推法,但精度也最差。
若一阶微分方程为
则
而
式中,h为步长。
显然步长越小,误差也小。
(2)改进的欧拉法是属于预测校正法,它将
左右两个点的斜率平均值作为
点处的计算。
即第一步计算预测值
第二步计算n+1点处的斜率
第三步计算校正值
(3)龙格—库塔法(Runge-KuttaMethod)一般形式为
均为待定常数,
下面介绍经典的三阶龙格—库塔法,表达式为:
下面介绍经典的四阶龙格—库塔法,表达式为:
%四阶M文件
function[x,y]=marunge4(dyfun,xspan,y0,h)
formatshort;
x=xspan
(1):
h:
xspan
(2);
y
(1)=y0;
forn=1:
(length(x)-1)
k1=feval(dyfun,x(n),y(n));
k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);
k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);
k4=feval(dyfun,x(n+1),y(n)+h*k3);
y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;
x=x'
y=y'
%三阶M文件
function[x,y]=marunge3(dyfun,xspan,y0,h)
k3=feval(dyfun,x(n)+h/2,y(n)+h/2*(2*k2-k1));
y(n+1)=y(n)+h*(k1+4*k2+k3)/6;
一阶微分方程
,初值条件为
,求
时的
向量。
在命令窗口输入如下程序:
dyfun=inline('
3*y/(1+x)'
[x,y]=marunge4(dyfun,[0,1],1,0.1);
[x'
y'
];
dsolve('
Dy=3*y/(1+x)'
'
y(0)=1'
x'
);
%用符号解,解上述微分方程式
x=[0:
0.1:
1];
y1=ans;
y1=1+3*x+3*x.^2+x.^3;
%微分方程式的精确解
err=y1-y'
%精确解与数值解之差。
plotyy(x,y,x,err)%绘制双纵坐标曲线,右侧坐标为y,左侧纵坐标为精确解与数值解之差
holdon
[x,y2]=marunge3(dyfun,[0,1],1,0.1);
y2'
err=y1-y2'
得到如图4.1所示:
图4.1微分方程dy/dx=3y/(1+x)精确解与数值解的误差曲线
实验五最小二乘法拟合多项式的实际应用
a.了解最小二乘拟合的基本原理和方法
b.掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法
c.通过实例学习如何用拟合方法解决实际问题,注意与插值方法的区别
d.了解各种参数辨识的原理和方法
e.通过范例展现由机理分析确定模型结构
设用一个m次多项式
----
(1)
来拟合一组数据
在节点
处,
的偏差设为
最小二乘法的基本思想是对所有数据点,拟合函数式
(1)的偏差
的平方和
-------
(2)
取最小值。
由于
为已知,故将
视作
的系数
的函数,不同的拟合多项式有不同的一组系数
,因而有不同的
值,即
于是上述数据拟合问题便归结为求多元函数的极值问题,欲使
取极小值,则
必须满足条件
对
(2)式求偏导数得
--------------(3)
由(3)等于0,有
---------------(4)
则(4)式可表示为
这是一个m+1阶对称的线性方程组,称为正规方程组,具有下面格式
此方程组的系数行列式不为零,故它有唯一解,将解得系数
代入
(1)式即可得所要求的拟合多项式。
functionli(x,y,m)
S=zeros(1,2*m+1);
T=zeros(m+1,1);
2*m+1
S(k)=sum(x.^(k-1));
m+1
T(k)=sum(x.^(k-1).*y);
A=zeros(m+1,m+1);
a=zeros(m+1,1);
forj=1:
A(i,j)=S(i+j-1);
a=A\T;
fprintf('
a[%d]=%f\n'
k,a(k));
本题数据选自2012年数学建模A题深圳人口预测,根据所给数据,预测深圳未来十年常住人口数量:
年份
年末常住人口(万人)
1979
31.41
1980
33.29
1981
36.69
1982
44.95
1983
59.52
1984
74.13
1985
88.15
1986
93.56
1987
105.44
1988
120.14
1989
141.6
1990
167.78
1991
226.76
1992
268.02
1993
335.97
1994
412.71
1995
449.15
1996
482.89
1997
527.75
1998
580.33
1999
632.56
在MATLAB中输入如下程序:
x=0:
1:
20;
y=[31.41,33.29,36.69,44.95,59.52,74.13,88.15,93.56,105.44,120.14,141.6,167.78,226.76,268.02,335.97,412.71,449.15,482.89,527.75,580.33,632.56];
li(x,y,3);
plot(x,y,'
*'
)
holdon;
x1=linspace(0,20);
y1=38.504224-4.881726*x1+1.787238*x1.^2;
plot(x1,y1,'
r'
li(x,y,5);
y2=46.220930-10.166090*x1+2.464142*x1.^2-0.022563*x1.^3;
plot(x1,y2,'
g'
li(x,y,9);
y3=21.266865+21.809792*x1-5.143542*x1.^2+0.578326*x1.^3-0.015022*x1.^4;
plot(x1,y3,'
b'
拟合曲线结果如图5.1所示:
图5.1多项式拟合
由上图可看出九次多项式拟合结果较好,则我们选用九次多项式进行预测;
总结
本文共进行了五项实验即列主元Gauss消去法、Newton插值、龙贝格求积公式、Runge-Kutta方法及最小二乘拟合的数值方法实际应用。
通过使用MATLAB数学软件,熟练的掌握了数学知识的运用,通过本次课设,充分理解了数值分析的某些作用以及意义所,并且通过本次课设,提高了我分析和解决问题的能力,做到学以致用。
参考文献
【1】杜延松,沈艳军等.数值分析及实验.北京:
科学出版社,2006.
【2】张德丰.Matlab数值分析与应用.北京:
国防工业出版社,2007.
【3】白峰衫.数值计算引论.北京:
高等教育出版社,2004.
【4】蔡大用.数值分析与实验学习指导.北京:
清华大学出版社,2001.
【5】关治,陆金甫.数值分析基础.北京:
高等教育出版社,1998.
【6】CleveB.Moler(著),喻文健(译).MATLAB数值计算.北京:
机械工业出版社,2006.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析 课程设计 数值 分析
![提示](https://static.bdocx.com/images/bang_tan.gif)