实验地点：研究生院机房 实验人：##

实验一 函数插值方法 

一、问题提出 

对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。

数据如下： 

（1） 

的值。（提示：结果为,   ） 

（2） 

二、要求 

1、 利用Lagrange插值公式 

编写出插值多项式程序； 

2、 给出插值多项式或分段三次插值多项式的表达式； 

3、 根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； 

4、 对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下： 

其中： 

三、目的和意义 

1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题； 

2、 明确插值多项式和分段插值多项式各自的优缺点； 

3、 熟悉插值方法的程序编制； 

4、 如果绘出插值函数的曲线，观察其光滑性。

四、实验学时：2学时

五、实验步骤： 

1．进入C或matlab开发环境；

2．根据实验内容和要求编写程序；

3．调试程序；

4．运行程序；

5．撰写报告,讨论分析实验结果.

六、实验程序及注释

1.程序一

function f=Lagrange(x,fx,inx) 

x=[0.4 0.55 0.65 0.8 0.95 1.05]

fx=[0.41075 0.57815 0.69675 0.90 1.0 1.25382]

inx=[0.596,0.99];

n=length(x);

m=length(inx); 

for i=1:m; 

   z=inx(i); 

   s=0.0; 

for k=1:n 

      p=1.0;  

for j=1:n 

if j~=k  

            p=p*(z-x(j))/(x(k)-x(j)); 

end 

      end  

s=p*fx(k)+s; 

end  

f(i)=s; 

end  

plot(x,fx,'O',inx,f)

运行结果：

x =

    0.4000    0.5500    0.6500    0.8000    0.9500    1.0500

fx =

    0.4108    0.5782    0.6967    0.9000    1.0000    1.2538

ans =

    0.6257    1.0542

图像：

2.程序二

function f=Lagrange(x,fx,inx) 

x=[1 2 3 4 5 6 7]

fx=[0.368 0.135 0.050 0.018 0.007 0.002 0.001]

inx=[1.8 6.15];

n=length(x);

m=length(inx); 

for i=1:m; 

   z=inx(i); 

   s=0.0; 

for k=1:n 

      p=1.0; 

for j=1:n 

if j~=k  

            p=p*(z-x(j))/(x(k)-x(j)); 

运行结果：

 x =

     1     2     3     4     5     6     7

fx =

    0.3680    0.1350    0.0500    0.0180    0.0070    0.0020    0.0010

ans =0.18  0.0013

图像：

实验二 函数逼近与曲线拟合 

一、问题提出 

从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。 

1、用最小二乘法进行曲线拟合； 

2、近似解析表达式为；

3、打印出拟合函数，并打印出与的误差，； 

4、另外选取一个近似表达式，尝试拟合效果的比较； 

5、* 绘制出曲线拟合图。 

三、目的和意义 

1、掌握曲线拟合的最小二乘法； 

2、最小二乘法亦可用于解超定线代数方程组； 

3、探索拟合函数的选择与拟合精度间的关系。

四、实验学时：2学时

五、实验步骤： 

1．进入C或matlab开发环境；

2．根据实验内容和要求编写程序；

3．调试程序；

4．运行程序；

5．撰写报告,讨论分析实验结果．

六、实验程序： 

x=[0 5 10 15 20 25 30 35 40 45 50 55];

y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.];

p=polyfit(x,y,3)

x1=0:2.5:55;

y1=polyval(p,x1);

plot(x,y,'*',x1,y1,'k')

七、实验结果： 

p =0.0000   -0.0052    0.2634    0.0178

七、实验图像：

实验三 数值积分与数值微分

一、基本题 

选用复合梯形公式，复合Simpson公式，Romberg算法，计算  

（1） 

（2） 

（3） 

   二、应用题

1.文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在两坐标轴上取天文测量单位（一天文单位为地球到太阳的平均距离：9300万里）在五个不同的时间对小行星作了五次观察，测得轨道上五个点的坐标数据如下表所示：         

现需要建立椭圆的方程以供研究。

 （1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组AX = b。

（2）用MATLAB求低价方程组的指令A / b求出待定系数 。

（3）卫星轨道是一个椭圆，其周长的计算公式是:

式中，a是椭圆的半长轴， 是地球中心与轨道中心（椭圆中心）的距离， 。其中h为近地点距离，H为远地点距离，R = 6371(km)为地球半径。

   有一颗人造卫星近地点距离h = 439 (km)，远地点距离H = 2384(km)。试分别按下列方案计算卫星轨道的周长，误差限取为 。

三、要求 

1、 编制数值积分算法的程序； 

2、 分别用两种算法计算同一个积分，并比较其结果； 

3、 分别取不同步长，试比较计算结果（如n = 10, 20等）； 

4、 给定精度要求ε，试用变步长算法，确定最佳步长。 

四、目的和意义 

1、 深刻认识数值积分法的意义； 

2、 明确数值积分精度与步长的关系； 

3、 根据定积分的计算方法，结合专业考虑给出一个二重积分的计算问题。

五、实验学时：2学时

六、实验步骤： 

1．进入C或matlab开发环境；

2．根据实验内容和要求编写程序；

3．调试程序；

4．运行程序；

5．撰写报告,讨论分析实验结果．

7、实验程序： 

clc

format long

clear

a=0;b=1;p=0;

n=4;m=0;

h=(b-a)/n;

for i=1:n-1

          x=0;

          x=a+h*i;

          m=m+sin(x)/x;

end

h=(b-a)/(2*n);

for k=1:2:2*n              

        p=p+sin(a+k*h)/(a+k*h);

end

Sn=h/3*[1+4*p+2*m+sin(b)/b] %辛普森公式结果

result=sinint(1)  %计算精确值

format short

八、实验结果：

 Sn =0.94608331088847

 result =0.94608307036718

实验四 常微分方程初值问题数值解法

一、基本题 

科学计算中经常遇到微分方程（组）初值问题，需要利用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. 小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2.

建立火箭升空过程的数学模型（微分方程）；

求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度．

2. 小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k，火箭升空过程的数学模型为

其中为火箭在时刻t的高度，m=1200-15t为火箭在时刻t的质量，T（=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t1 (=900/15=60秒)为引擎关闭时刻．

今测得一组数据如下（t~时间（秒），x ~高度（米），v~速度（米/秒））：

   1．用每一个数据(t,x,v)计算一个k的估计值（共11个），再用它们来估计k。

   2．用这组数据拟合一个k　．

请你分别用这两种方法给出k的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）．

三、要求 

1、 根据初值问题数值算法，分别选择二个初值问题编程计算； 

2、 试分别取不同步长，考察某节点处数值解的误差变化情况； 

3、 试用不同算法求解某初值问题，结果有何异常； 

4、 分析各个算法的优缺点。 

四、目的和意义 

1、 熟悉各种初值问题的算法，编出算法程序； 

2、 明确各种算法的精度与所选步长有密切关系； 

3、 通过计算更加了解各种算法的优越性。 

五、实验学时：2学时

六、实验步骤： 

1．进入C或matlab开发环境；

2．根据实验内容和要求编写程序；

3．调试程序；

4．运行程序；

5．撰写报告,讨论分析实验结果.

    七、实验程序：

#include

#include 

#define f(x,y)(y-2*x/y) 

int main() 

{   int m;  int i;

    double a,b,y0; 

    double xn,yn,yn1; 

    double k1,k2,k3,k4; double h;

    printf("\\nInput the begin and end of x:"); 

    scanf("%lf%lf",&a,&b);  printf("nInput the y value at %f:",a);

scanf("%lf",&y0);  printf("nInput m value[divide(%f,%f)]:",a,b);

scanf("%d",&m); if(m<=0) {  printf("Please input a number larger than 1.\\n"); return 1; }

h=(b-a)/m; xn=a;yn=y0; 

for(i=1;i<=m;i++) {  k1=f(xn,yn); 

    k2=f((xn+h/2),(yn+h*k1/2));

    k3=f((xn+h/2),(yn+h*k2/2)); 

    k4=f((xn+h),(yn+h*k3)); 

    yn1=yn+h/6*(k1+2*k2+2*k3+k4); 

    xn+=h;

    printf("x%d=%f,y%d=%f\\n",i,xn,i,yn1); yn=yn1; } 

    scanf("%lf",yn);}下载本文