数值分析实验报告.docx
- 文档编号:23979792
- 上传时间:2023-05-23
- 格式:DOCX
- 页数:19
- 大小:215.86KB
数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告.docx(19页珍藏版)》请在冰豆网上搜索。
数值分析实验报告
机电工程学院机械工程陈星星6720150109
《数值分析》课程设计实验报告
实验一函数插值方法
一、问题提出
对于给定的一元函数
的n+1个节点值
。
试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。
数据如下:
(1)
0.4
0.55
0.65
0.80
0.95
1.05
0.41075
0.57815
0.69675
0.90
1.00
1.25382
求五次Lagrange多项式
,计算
的值。
(提示:
结果为
)
实验步骤:
第一步:
先在matlab中定义lagran的M文件为拉格朗日函数
代码为:
function[c,l]=lagran(x,y)
w=length(x);
n=w-1;
l=zeros(w,w);
fork=1:
n+1
v=1;
forj=1:
n+1
if(k~=j)
v=conv(v,poly(x(j)))/(x(k)-x(j));
end
end
l(k,:
)=v;
end
c=y*l;
end
第二步:
然后在matlab命令窗口输入:
>>>>x=[0.40.550.650.80,0.951.05];y=[0.410750.578150.696750.901.001.25382];
>>p=lagran(x,y)
回车得到:
P=
121.6264-422.7503572.5667-377.2549121.9718-15.0845
由此得出所求拉格朗日多项式为
p(x)=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845
第三步:
在编辑窗口输入如下命令:
>>x=[0.40.550.650.80,0.951.05];
>>y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845;
>>plot(x,y)
命令执行后得到如下图所示图形,然后
>>x=0.596;
>>y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084
y=0.6257得到f(0.596)=0.6257同理得到f(0.99)=1.0542
(2)
1
2
3
4
5
6
7
0.368
0.135
0.050
0.018
0.007
0.002
0.001
试构造Lagrange多项式
,和分段三次插值多项式,计算的
值。
(提示:
结果为
)
实验步骤:
第一步定义
function[c,l]=lagran(x,y)
w=length(x);
n=w-1;
l=zeros(w,w);
fork=1:
n+1
v=1;
forj=1:
n+1
if(k~=j)
v=conv(v,poly(x(j)))/(x(k)-x(j));
end
end
l(k,:
)=v;
end
c=y*l;
end
定义完拉格朗日M文件
第二步:
然后在matlab命令窗口输入:
>>x=[1234567];y=[0.3680.1350.0500.0180.0070.0020.001];
>>p=lagran(x,y)
回车得到:
由此得出所求拉格朗日多项式为
p(x)=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950
第三步:
在编辑窗口输入如下命令:
>>x=[1234567];
>>y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950;
>>plot(x,y)
命令执行后得到如下图所示图形,然后
>>x=1.8;
>>y=4304240283865561*x.^6/73786976294838206464-7417128346304051*x.^5/4611686018427387904+223*x.^4/12000-2821*x.^3/24000+994976512675275*x.^2/2251799813685248-19367*x./20000+199/200
y=
0.1648得到f(1.8)=0.1648同理得到f(6.15)=0.0013
实验结论:
插值是在离散数据的基础上补插连续函数,使得这条连续曲线通过全部给定的离散数据点,它是离散函数逼近的重要方法,利用它可通过函数在有限个点处的取值状况,估算出函数在其他点处的近似值。
实验二函数逼近与曲线拟合
一、问题提出
从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量与时间t的拟合曲线。
t(分)
0510152025303540455055
01.272.162.863.443.874.154.374.514.584.024.64
第一步先写出线性最小二乘法的M文件
functionc=lspoly(x,y,m)
n=length(x);
b=zeros(1:
m+1);
f=zeros(n,m+1);
fork=1:
m+1
f(:
k)=x.^(k-1);
end
a=f'*f;
b=f'*y';
c=a\b;
c=flipud(c);
第二步在命令窗口输入:
>>lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)
回车得到:
ans=
-0.0024
0.2037
0.2305
即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305
在编辑窗口输入如下命令:
>>x=[0,5,10,15,20,25,30,35,40,45,50,55];
>>y=-0.0024*x.^2+0.2037*x+0.2305;
>>plot(x,y)
命令执行得到如下图
实验结论:
分析复杂实验数据时,常采用分段曲线拟合方法。
利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。
分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。
实验四线方程组的直接解法
一、问题提出
给出下列几个不同类型的线性方程组,请用适当算法计算其解。
1、设线性方程组
2、设对称正定阵系数阵线方程组
3、三对角形线性方程组
实验步骤:
列主元高斯消去法的matlab的M文件程序
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
error('ThecolumnsofAmustbeequalthelengthofb!
');
return;
end
%开始计算,先赋初值
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
end
ifa_max<1e-10
index=0;return;
end
%交换两行
ifr>k
forj=k:
n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;det=-det;
end
%消元过程
fori=k+1:
n
m=A(i,k)/A(k,k);
forj=k+1:
n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<1e-10
index=0;return;
end
fork=n:
-1:
1
forj=k+1:
n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
然后在命令窗口输入
>>A=[42-3-1210000;86-5-3650100;42-2-132-1031;0-215-13-1194;-426-167-3323;86-8571726-35;02-13-425301;1610-11-917342-122;462-713920124;00-18-3-24-863-1];
>>b=[51232346133819-21];
>>gauss(A,b)
ans=
1.0000
-1.0000
0.0000
1.0000
2.0000
0.0000
3.0000
1.0000
-1.0000
2.0000
高斯-约当消去法maltab的M文件程序
function[x,flag]=Gau_Jor(A,b)
%求线形方程组的列主元Gauss-约当法消去法,其中,
%A为方程组的系数矩阵;
%b为方程组的右端项;
%x为方程组的解;
[n,m]=size(A);nb=length(b);
%当方程组行与列的维数不相等时,停止计算,并输出出错信息。
ifn~=m
error('TherowsandcolumnsofmatrixAmustbeequal!
');
return;
end
%当方程组与右端项的维数不匹配时,停止计算,并输出出错信息
ifm~=nb
error('ThecolumnsofAmustbeequalthelengthofb!
');
return;
end
%开始计算,先赋初值
flag='ok';x=zeros(n,1);
fork=1:
n
%选主元
max1=0;
fori=k:
n
ifabs(A(i,k))>max1
max1=abs(A(i,k));r=i;
end
end
ifmax1<1e-10
falg='failure';return;
end
%交换两行
ifr>k
forj=k:
n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;
end
%消元过程
b(k)=b(k)/A(k,k);
forj=k+1:
n
A(k,j)=A(k,j)/A(k,k);
end
fori=1:
n
ifi~=k
forj=k+1:
n
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
b(i)=b(i)-A(i,k)*b(k);
end
end
end
%输出x
fori=1:
n
x(i)=b(i);
end
然后保存后在命令窗口输入:
>>A=[42-402400;22-1-21320;-4-1141-8-356;0-216-1-4-33;21-8-1224-10-3;43-3-44111-4;025-3-101142;0063-3-4219];
>>b=[0-620239-22-1545];
>>Gau_Jor(A,b)
ans=
121.1481
-140.1127
29.7515
-60.1528
10.9120
-26.7963
5.4259
-2.0185
实验结论:
用LU法,调用matlab中的函数lu中,L往往不是一个下三角,但可以直接计算不用它的结果来计算,不用进行行变换。
如果进行行变b也要变,这样会很麻烦。
实验八常微分方程初值问题数值解法
一、基本题
科学计算中经常遇到微分方程(组)初值问题,需要利用Euler法,改进Euler法,Rung-Kutta方法求其数值解,诸如以下问题:
(1)
分别取h=0.1,0.2,0.4时数值解。
初值问题的精确解
。
(2)
用r=3的Adams显式和预-校式求解
取步长h=0.1,用四阶标准R-K方法求值。
(3)
用改进Euler法或四阶标准R-K方法求解
取步长0.01,计算
数值解,参考结果
。
(4)利用四阶标准R-K方法求二阶方程初值问题的数值解
(I)
(II)
(III)
(IV)
实验步骤:
function [x,y]=euler(fun,x0,xfinal,y0,n);
ifnargin<5,n=50;
end
h=(xfinal-x0)/n;
x
(1)=x0;y
(1)=y0;
fori=1:
n;
x(i+1)=x(i)+h;
y(i+1)=y(i)+h*feval(fun,x(i),y(i));
end
实验程序及分析
(Ⅰ)
(1)、算法程序
functionE=Euler_1(fun,x0,y0,xN,N)
%Euler向前公式,其中
%fun为一阶微分方程的函数
%x0,y0为初始条件
%xN为取值范围的一个端点
%h为区间步长
%N为区间个数
%x为Xn构成的向量
%y为yn构成的向量
x=zeros(1,N+1);y=zeros(1,N+1);
x
(1)=x0;y
(1)=y0;
h=(xN-x0)/N;
forn=1:
N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));
end
T=[x',y']
functionz=f(x,y)
z=4*x/y-x*y;
(2)、运行程序
>>Euler_1('f',0,3,2,20)
四、实验结果
>>Euler_1('f',0,3,2,20)
T=
03.0000
0.10002.9836
0.20002.9517
0.30002.9058
0.40002.8481
0.50002.7810
0.60002.7073
0.70002.6297
0.80002.5511
0.90002.4739
1.00002.4004
1.10002.3325
1.20002.2714
1.30002.2177
1.40002.1717
1.50002.1332
1.60002.1017
1.70002.0765
1.80002.0567
1.90002.0414
2.00002.0299
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告