江西理工大学研究生学院
《数值分析》课程设计实验报告
姓 名： 彭剑
学 号： 6120130166
专 业： 结构工程
指导老师： 岳老师
实验一 函数插值方法
一、问题提出
二、实验程序及注释
第一步：现在matlab中定义lagran的M文件为拉格朗日函数代码为：
function[c,l]=lagran(x,y)
w=length(x);
n=w-1;
l=zeros(w,w);
for k=1:n+1
v=1;
for j=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.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382];
>>lagran(x,y)
回车得到：
ans =
121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845
由此得出所求拉格朗日多项式为
p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845
第三步：在编辑窗口输入如下命令：
>> x=[0.4 0.55 0.65 0.80,0.95 1.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.6262
得到f（0.596）=0.6262 同理得到f（0.99）=1.0547
三、实验数据结果
同理：
第一步：定义拉格朗日M文件同上题（上题已定义这里直接引用）
第二步：然后在matlab命令窗口输入：
>>>> x=[1 2 3 4 5 6 7]; y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001];
>> ans=lagran(x,y)
回车得到：
ans =
0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950
由此得出所求拉格朗日多项式为
p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950
第三步：在编辑窗口输入如下命令：
>> x=[1 2 3 4 5 6 7];
>> 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=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;
y =0.1670
得到f（1.8）=0.6270 同理得到f（6.15）=2.3644。
四、实验结论：
插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值
实验二 函数逼近与曲线拟合
二、实验程序及注释
第一步先写出线性最小二乘法的M文件
function c=lspoly(x,y,m)
n=length(x);
b=zeros(1:m+1);
f=zeros(n,m+1);
for k=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);
% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。
if n~=m
error('The rows and columns of matrix A must be equal!');
return;
end
% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息
if m~=nb
error('The columns of A must be equal the length of b!');
return;
end
% 开始计算，先赋初值
index=1;det=1;x=zeros(n,1);
for k=1:n-1
% 选主元
a_max=0;
for i=k:n
if abs(A(i,k))>a_max
a_max=abs(A(i,k));r=i;
end
end
if a_max<1e-10
index=0;return;
end
% 交换两行
if r>k
for j=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
% 消元过程
for i=k+1:n
m=A(i,k)/A(k,k);
for j=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);
% 回代过程
if abs(A(n,n))<1e-10
index=0;return;
end
for k=n:-1:1
for j=k+1:n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
然后在命令窗口输入
>> A=[4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1];
>> b=[5 12 3 2 3 46 13 38 19 -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);
% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。
if n~=m
error('The rows and columns of matrix A must be equal!');
return;
end
% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息
if m~=nb
error('The columns of A must be equal the length of b!');
return;
end
% 开始计算，先赋初值
flag='ok';x=zeros(n,1);
for k=1:n
% 选主元
max1=0;
for i=k:n
if abs(A(i,k))>max1
max1=abs(A(i,k));r=i;
end
end
if max1<1e-10
falg='failure';return;
end
% 交换两行
if r>k
for j=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);
for j=k+1:n
A(k,j)=A(k,j)/A(k,k);
end
for i=1:n
if i~=k
for j=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
for i=1:n
x(i)=b(i);
end
然后保存后在命令窗口输入：
>> A=[4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1];
>> b=[5 12 3 2 3 46 13 38 19 -21];
>> Gau_Jor(A,b)
ans =
1.0000
-1.0000
0.0000
1.0000
2.0000
0.0000
3.0000
1.0000
-1.0000
2.0000
在利用两种方法都互相检验，验证是结果误差很小，其他几题同理可得。
三、实验结论：
用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) ?
二、要求
1、 根据初值问题数值算法，分别选择二个初值问题编程计算；
2、 试分别取不同步长，考察某节点处
数值解的误差变化情况；
3、 试用不同算法求解某初值问题，结果有何异常；
4、 分析各个算法的优缺点。
三、实验程序及分析
（Ⅰ）
（1）、算法程序
function E =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;
for n=1:N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));
end
T=[x',y']
function z=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 =
0 3.0000
0.1000 2.9836
0.2000 2.9517
0.3000 2.9058
0.4000 2.8481
0.5000 2.7810
0.6000 2.7073
0.7000 2.6297
0.8000 2.5511
0.9000 2.4739
1.0000 2.4004
1.1000 2.3325
1.2000 2.2714
1.3000 2.2177
1.4000 2.1717
1.5000 2.1332
1.6000 2.1017
1.7000 2.0765
1.8000 2.0567
1.9000 2.0414
2.0000 2.0299
五、实验结论：
欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较低。改进欧拉格式较欧拉格式提高了精度，其截断误差比欧拉格式提高了一阶。龙格－库塔法是用于模拟长微分方程的解的重要的一类隐式或显式迭代法。工程中很多的地方用到龙格库塔求解微分方程的数值解，龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。