Matlab for浙江大学1-6概论

简单实用,如数学,物理中的特殊函数表,数 理统计中用的概率分布表等.
三次样条插值
n
概念 s ( x ) 满足: a) 在每个小区间 [ xi1 , xi ](i = 1,..., n) 上是3 次多项式插值; b) 在 a ≤ x ≤ b 上二阶导数连续; c) s ( xi ) = yi , i = 0,1,..., n
y=Ex_3(t)
a=8755;b=6810; y=sqrt(a^2*sin(t).^2+b^2*cos(t).^2); 在Matlab工作区键入以下程序: t=0:pi/10:pi/2; y1=Ex_3(t); L1=4*trapz(t,y1); L2=4*quad('Ex_3',0,pi/2,1e-6) 用长格式输出结果为: L1=4.908996526785276e+004 L2=4.908996531830460e+004
矩形公式: Ln = h∑ f k
k =0 n
n1
Rn = h∑ f k
k =1
n
Matlab上的实现 输入: a b n y h=(b-a)/n L=sum(y(1:n))*h R=sum(y(2:n+1))*h Z=cumsum(y) L=z(n)*h, R=(z(n+1)-z(1))*h
梯形求积公式
…… ……
xn yn
拉格朗日插值
n
Matlab上的实现 y 输入值:n 个节点对应数组 x0 , 0,以及 m 个插值点对应的数组 x 输出值:m 个插值对应的数组 y lagr.m y=lagr(x0,y0,x)
y=lagr(x0,y0,x)
m=length(X0); N=zeros(m,1); y=0; for i=1:m N(i)=1; for j=1:m if j~=i N(i)=N(i)*(x-X0(j))/(X0(i)-X0(j)); end end y=y+y0(i)*N(i); end
n
Matlab上的实现 y 输入值: n 个节点对应数组 x0 , 0,以及 m 个插值点对应的数组 x 输出值:m 个插值对应的数组 y y=interp1(x0,y0,x,'spline') y=spline(x0,y0,x) Runge现象的三次样条插值如前图
n
三次样条插值总结 优点:具有良好的收敛性,节点处的光滑 性提高; 缺点:误差估计难.
y
S=quad('fun',a,b)
tol=1e-3
S=quad('fun',a,b,tol)
n
误差估计 f ∈ C 4 ( a, b) ,M 4 = max f (4) ( x) 则有 若
h4 R( f , Sn ) = I Sn ≤ M 4 (b a) 180 = O(h4 )
四阶收敛,代数精度为3,精度高,但 计算工作量高.
二,数值积分
n
被积函数 f ( x ) 的原函数 F ( x ) 无法由基本初 等函数经过四则及复合运算构成.如
n
sin x ∫ e dx , ∫ x dx f ( x )是由测量或数值计算给出的一张表.
x2
x y
x0 y1
x1 y2
x2 y3
…… ……
xn yn
矩形求积公式
n
定义
a = x0 < x1 < L < xk < L < xn = b ba h= , f k = f ( xk ), k = 0,1,K, n n
Ex_1.m
x0=[0 3 5 7 9 11 12 13 14 15]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6]; x=0:0.1:15; y1=lagr(x0,y0,x); y2=interp1(x0,y0,x); y3=interp1(x0,y0,x,'spli ne'); [x',y1',y2',y3'] subplot(3,1,1) plot(x0,y0,'+k',x,y1,'r') grid title('lagrange') subplot(3,1,2) plot(x0,y0,'+k',x,y2,'r') grid title('piecewise linear') subplot(3,1,3) plot(x0,y0,'+k',x,y3,'r') grid title('spline')
三,常微分方程数值解
n
问题的提法:设有一阶方程和初始条件 y′ = f ( x, y ), y ( x0 ) = y0 其中 f 适当光滑,对 y 满足Lipschitz条件, 以保证方程的解存在且唯一. 该微分方程的解析解 y = y ( x ) 难以求解或根本无 解析解,而是在一系列离散点 x0 < x1 < L < xn 上 求 y ( xn ) 的近似值 yn ( n = 0,1,L, n) 通常取等步长 h ,即 xn = x0 + nh( n = 1, 2,L)
n
Matlab上的实现 y 输入值: n 个节点对应数组 x0 , 0,以及 m 个插值点对应的数组 x 输出值:m 个插值对应的数组 y y=interp1(x0,y0,x) Runge现象的分段线性插值如前图
n
分段线性插值总结 优点:简单实用,计算 x 点的插值时,只用 到 x 左右的两个节点,计算量与节 点个数 n 无关,n 越大,分段越多 插值误差越小.保证收敛性. 缺点:在节点处不光滑,精度差些.
向前欧拉公式
n
定义 y(x ) y(x ) 在小区间 [ xn , xn+1 ] 上用差商 代替方程左 h 端的导数 y′ ,而方程右端已知函数 f ( x, y )中的 取左端点 xn,可得 y ( xn+1 ) ≈ y ( xn ) + hf ( xn , y ( xn )) 从 x0 出发,由初值 y ( x0 ) = y0 代入得到 y ( x1 )的近似 值 y1 ,继续代入右端,得到向前欧拉公式如下: yn +1 = yn + hf ( xn , yn ), n = 0,1,L
卫星轨道长度计算范例
n
人Байду номын сангаас地球卫星轨道可视为平面上的椭圆, 我国第一颗人造地球卫星近地点距地球表 面439cm,远地点距地球表面2384cm, 地球半径为6371cm,求该卫星的轨道长 度.
Ex_3.m
解题思路
n
卫星轨道椭圆的参数方程为
x = a cos t , y = b sin t (0 ≤ t ≤ 2π )
n
分段线性插值公式
x x j 1 x x j 1 j x x j +1 l j ( x) = x j x j +1 0 x j 1 ≤ x ≤ x j , j ≠ 0 x j ≤ x ≤ x j +1 , j ≠ n [ a, b] [ x j 1 , x j +1 ]
n
拉格朗日多项式插值方法总结 优点:高次多项式插值,插值曲线光滑, 误差估计有表达式; 缺点:有振荡现象,收敛性不能保证.
用于理论分析,实际意义不大
分段线性插值
n
插值多项式的振荡 n ↑ ,Rn ( x) ↓ 同时 Ln ( x) 的光滑性变坏,有 时会出现很大的振荡. 龙格(Runge)现象:
n +1 n
隐式公式,不能直接计算 yn+1 精度为1阶
梯形公式
n
为了提高精度到2阶,将向前和向后欧拉公式 加以平均,得到梯形公式: h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , yn +1 )], n = 0,1,L 2 隐式公式,不能直接计算 yn+1 精度为2阶
n
定义
a = x0 < x1 < L < xk < L < xn = b ba h= , f k = f ( xk ), k = 0,1,K, n n n1 h 梯形公式: Tn = h∑ f k + ( f 0 + f n ) 2 k =1
n
Matlab上的实现 输入: a b n y h=(b-a)/n T=trapz(y)*h(输入数组x,输出为积 分值,但步长为单位步长) T=trapz(x,y)*h(输入x,y为同长度的 数组,输出y对x的积分,但步长不一定 相等)
简单实用,应用广泛.
机床加工范例
n
待加工零件的外形根据工艺要求由一组数 据 ( x, y ) 给出,需要得到 x 坐标每改变 0.1 时的 y 坐标.试完成加工所需数据,画出 曲线,并求出 x = 0 处的曲线斜率和 13 ≤ x ≤ 15 范围内 y 的最小值. Ex_1.m
x 0 3 5 7 9 11 12 13 14 15 y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6
确定 建立 数值 实际 → 数学 → 计算 问题 模型 方法
编制 程序 → 上机 计算 结果
一,插值
n
a = x0 < x1 < ... < xn = b 是 [ a, b] 上 n + 1 个互异 x* ( ≠ x j ) 处的插值 y* 的点,求任一插值点 x y x0 y1 x1 y2 x2 y3
定义
a = x0 < x1 < L < xk < L < xn = b, ba h= , n f k = f ( xk ), m = n/2 k = 0,1,K, n
m1 m 1 h Sn = ( f 0 + f n + 4∑ f 2 k +1 + 2∑ f 2 k ) 3 k =0 k =1
n
Matlab上的实现 输入: a b n = 2m h=(b-a)/n fun.m
2π
a, b 分别是长,短半轴.根据计算参数方程的
弧长的公式,椭圆长度可表为如下积分: 1
L = 4∫ ( a sin t + b cos t ) 2 dt
2 2 2 2 0
称为椭圆积分,它无法用解析方法计算,根 据所给数据a=6371+2384=8755, b=6371+439=6810.下面是用梯形公式和辛 普森公式计算的程序:
y=eulerpro(x,y0)
n=length(x); h=(x(n)-x(1))/(n-1); y(1)=y0; for i=1:n-1 y1(i+1)=y(i)+h*(y(i)+2*x(i)); y(i+1)=y(i)+h/2*((y(i)+2*x(i)+y1(i +1)+2*x(i+1))); end
n
误差估计 f ∈ C 2 ( a, b) , M 2 = max f ′′( x ) 则有 若
h2 R ( f , Tn ) = I Tn ≤ M 2 (b a ) 12 = O(h2 )
二阶收敛,代数精度为1,精度不高, 但计算工作量低,简单,运行方便.
辛普森(Simpson)求积公式
n
改进欧拉公式
n
定义 在小区间 [ xn , xn+1 ] 上先由向前欧拉公式算出 yn+1 的预测值 yn+1 ;再把它代入梯形公式右端,作 为校正,即 yn +1 = yn + hf ( xn , yn ) h yn +1 = yn + [ f ( xn , yn ) + f ( xn+1 , yn +1 )], n = 0,1,L 2 精度为2阶 Eulerpro.m
1 g ( x) = 2 1+ x 5 ≤ x ≤ 5
Runge.m
Runge(n)
x0=[-5:1:5]; y0=1./(1+x0.^2); x=[-5:0.1:5]; y=lagr(x0,y0,x); y1=1./(1+x.^2); y2=interp1(x0,y0,x); y3=interp1(x0,y0,x,'spli ne'); plot(x,y,'--r',x,y1,'g',x,y2,'-m',x,y3,'-b')
n +1 n
精度为1阶 Euler.m
y=Euler(x,y0)
y(1)=y0; n=length(x); h=(x(n)-x(1))/(n-1); for i=1:n-1 y(i+1)=y(i)+h*(y(i)+2*x(i)); end
向后欧拉公式
n
定义 y(x ) y(x ) 在小区间 [ xn , xn+1 ] 上用差商 代替方程左 h 端的导数 y′ ,而方程右端已知函数 f ( x, y )中的 取右端点xn+1 ,可得 y ( xn+1 ) ≈ y ( xn ) + hf ( xn+1 , y ( xn+1 )) 从 x0 出发,由初值 y ( x0 ) = y0 代入得到 y ( x1 )的近似 值 y1 ,继续代入右端,得到向后欧拉公式如下: yn +1 = yn + hf ( xn +1 , yn +1 ), n = 0,1,L