2016/3/18 8
• 3-37
• • • • • • • • • • • • • x = 0:10; % 给出已知基准数据 y = sin(x); xi = 0:.25:10; % 在两个基准数据点间插入3个点 yi1 = interp1(x,y,xi,'*nearest'); yi2 = interp1(x,y,xi,'*linear'); yi3 = interp1(x,y,xi,'*spline'); yi4 = interp1(x,y,xi,'*cubic'); % 分别用四种方法求中间插值 plot(x,y,'o',xi,yi1, 'r:',xi,yi2, 'g-',xi,yi3,'k.-',xi,yi4, 'm--') % 绘图，并标注各曲线代表的插值方法，如图3-8 legend('原始数据','最近点插值','线性插值','样条插值','立方插值')
12]内的数值解，且满足初始条件
y1 0 0 y2 0 1 y 0 1 3
15
2016/3/18
• 例3-39
• • • • • • • • • • • % 首先编写方程函数，名为rigid.m function dy = rigid(t,y) dy = zeros(3,1); % 生成3×1的矩阵 dy(1) = y(2) * y(3); % 第一个元素是 dy(2) = -y(1) * y(3); % 第二个元素是 dy(3) = -0.51 * y(1) * y(2); % 第三个元素是
例：根据经验公式y=a+bx2，拟合如下 数据：
xi 19 25 31 38 44
yi
19.0
32.3
49.0
73.3
98.8
x=[19 25 31 38 44]; y=[19.0 32.3 49.0 73.3 98.8]; x1=x.^2; x1=[ones(5,1),x1']; ab=x1\y' ab =
4
2016/3/18
例：拟合以下数据
x y 0.5 1.75 1.0 2.45 1.5 3.81 2.0 4.80 2.5 8.00 3.0 8.60
x=[0.5 1.0 1.5 2.0 2.5 3.0]; y=[1.75 2.45 3.81 4.80 8.00 8.60]; a=polyfit(x,y,2)
3
2016/3/18
• 例3-36
• • • • • • • • • • • • • • • x = (0: 0.1: 2.5)'; % 给出一组数据，为误差函数的一个区间 y = erf(x); p = polyfit(x,y,6) % 计算该区间内6阶拟合多项式的向量
x1 = (0: 0.1: 5)'; % 将区间增长一倍 y1 = erf(x1); % 计算误差函数在新区间内的函数值 f = polyval(p,x1); % 计算6阶拟合曲线在新区间内的取值 plot(x1,y1,'o',x1,f,'-') % 绘图，比较前后区间内的曲线拟合效果，如图3-7 axis([0 5 0 2])
x1= 0:0.01:1; f2 = polyval(y2,x1); % 2阶拟合曲线在各点的函数值 y10 = polyfit(x,y,10) % 计算10阶拟合的多项式向量 f10 = polyval(y10,x1); % 10阶拟合曲线在各点的函数值 plot(x,y,'o',x1,f2,':',x1,f10,'k')
y2 y1 y3 y2 原方程变为： y y 1 2 2 y 2 y 4 1 3 2 3
y1 0 1 初始条件： y2 0 1 y 0 0 3
2016/3/18
19
• 例3-40
• 训练任务：请编制名为odet3.m方程函数 • tspan = [0 12]; • % 在同一目录下，计算方程数值解。给出时间区间和初始 值 • y0 = [-1 1 0]; • [t,Y] = ode45('odet3',tspan,y0); • % 采用ode45算法 • plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.') • % 绘制计算结果并标注，如图3-15 • legend('Y1','Y2','Y3')
2016/3/18
2
• 例3-35
• • • • • • • • • • • • • • • x = (0:0.1:1)'; y = [-0.4 1.9 3.2 6.2 7.1 7.3 7.7 9.6 9.5 9.3 12]'; % 给出一组11个点数据 y2 = polyfit(x,y,2) % 计算2阶拟合的多项式向量
ZI4 = interp2(X,Y,Z,XI,YI,'*cubic'); % 计算用cubic法所得二维插值 figure(5) mesh(XI,YI,ZI4) % 绘制双立方插值法曲面图，如图3-13
13
2016/3/18
3.9常微分方程初值问题的数值解法 龙格－库塔法简介 龙格－库塔法的实现 基于龙格－库塔法，MATLAB提供了求常微分方程数值解 的函数，一般调用格式为： [t,y]=ode23('fname',tspan,y0) [t,y]=ode45('fname',tspan,y0) 其中fname是定义f(t,y)的函数文件名，该函数文件必须返回 一个列向量。tspan形式为[t0,tf],表示求解区间。y0是初始 状态列向量。t和y分别给出时间向量和相应的状态向量。
2016/3/18 14
• 一阶常微分方程的初值问题求解
y f x, y , a x b y x0 y 0
y 2 y3 y1 •例3-39，求微分方程组 y1 y3 在区间[0 y2 y 0.51y y 1 2 3
ZI2 = interp2(X,Y,Z,XI,YI,'*linear'); % 计算用linear法所得二维插值 figure(3) mesh(XI,YI,ZI2) % 绘制双线性插值法曲面图，如图3-11
12
2016/3/18
ZI3 = interp2(X,Y,Z,XI,YI,'*spline'); • % 计算用spline法所得二维插值 • figure(4) • mesh(XI,YI,ZI3) • % 绘制双样条插值法曲面图，如图3-12 • • • • •
2016/3/18
21
例设有初值问题，试求其数值解，并与精确解相比较 (精确解为y(t)=)。
(1) 建立函数文件funt.m。 function yp=funt(t,y) yp=(y^2-t-2)/4/(t+1); (2) 求解微分方程。 t0=0;tf=10; y0=2; [t,y]=ode23('funt',[t0,tf],y0); %求数值解 y1=sqrt(t+1)+1; %求精确解 t' y' y1' 2016/3/18 y为数值解，y1为精确值，显然两者近似。
10 8
a =
0.4900 1.2501 0.8560
6 4 2 0 0.5
x1=[0.5:0.05:3.0]; y1=a(3)+a(2)*x1+a(1)*x1.^2; plot(x,y,'*') hold on plot(x1,y1,'-r')
2016/3/18
1
1.5
2
2.5
3
5
• Matlab提供了一维、二维、 三次样条等
许多插值选择
2016/3/18 7
• 1.一维插值 yi=interp1(x,Y,xi,method)：求已知同维数据x和Y， 运用method指定的方法计算插值点xi处的数值yi。 method有如下四种： nearest 最近点差值，取与最近的已知数据点的值 linear 线性差值，直线连接数据点，插值点值位于 直线上 spline 样条插值，用三次样条曲线通过数据点，根 据曲线进行差值 cubic 立方差值，用三次曲线拟合并通过数据点
2016/3/18 17
• 高阶微分方程的初值问题求解
F y, y, y,, yn , x 0
高阶微分方程：


n 1 , yn 1 0 y0 y0 y0 , y0 y0
把高阶微分方程转换为一阶方程组，初始条件也做相应的替换。通常令： n 1 上式改写为： 1 2 n
9
2016/3/18
2.二维插值 ZI=interp2(X,Y,Z,XI,YI,method)：已知同维数 据X,Y和Z，运用method指定的方法计算插 值点(XI,YI)处的数值ZI。
2016/3/18
10
• 例3-38
• • • • • • • • • • • axis([-3 3 -3 3 -5 20]) % 确定坐标轴的范围 [X,Y] = meshgrid(-3:.5:3); % 生成网格 Z = peaks(X,Y); % 计算peaks函数值 figure(1) mesh(X,Y,Z) % 绘制原始数据曲面图，如图3-9
100
80
60