1 / 59

Chapter 7 MATLAB 在计算方法中的应用 7.1 插值 7.1.1 Lagrange 插值 对给定的 n 个自变量 x1, x2, …, xn 及对应的

Chapter 7 MATLAB 在计算方法中的应用 7.1 插值 7.1.1 Lagrange 插值 对给定的 n 个自变量 x1, x2, …, xn 及对应的 函数值 y1, y2, …, yn, 求插值区间内任意 x 的函数 值 y。 lagrange 插值公式： y(x)= ∑ y k （ ∏ ）. x-x j. n. n. x k -x j. j=1. k=1. j≠k. 将上述插值公式编写成 M 文件，取名 lagrange.m 句法结构：

oren-ochoa
Download Presentation

Chapter 7 MATLAB 在计算方法中的应用 7.1 插值 7.1.1 Lagrange 插值 对给定的 n 个自变量 x1, x2, …, xn 及对应的

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Chapter 7 MATLAB在计算方法中的应用 7.1 插值 7.1.1 Lagrange 插值 对给定的n个自变量x1, x2, …, xn及对应的 函数值y1, y2, …, yn, 求插值区间内任意x的函数 值y。 lagrange 插值公式： y(x)= ∑ yk（∏ ） x-xj n n xk-xj j=1 k=1 j≠k

  2. 将上述插值公式编写成M文件，取名lagrange.m 句法结构： function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n

  3. if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end

  4. 例1：x=0:0.2:2*pi; y=sin(x);x1=0:0.1:2*pi; 插值： y1=lagrange(x,y,x1); 精确值： y2=sin(x1); 检查插值效果：plot(x,y,’r’,x1,y1,’bo’) plot(x1,y1-sin(x1)) 例2：x=[-5:1:5]; y=1./(1+x.^2);x0=[-5:0.1:5]; 插值： y0=lagrange(x,y,x0); 精确值： y1= 1./(1+x0.^2); 检查插值效果：plot(x0,y0,’--r’,x0,y1,’-b’)

  5. 7.1.2 一维插值命令 interp1 句法结构：y=interp1(x0,y0,x) 例：将上例中的lagrange换成interp1,其他相同。 nearest:线性最近项插值，linear:线性插值 spline:三次样条插值，cubic：三次插值。

  6. 例 正弦曲线的插值示例。 >>x=0:0.1:10; >>y=sin(x); >>xi=0:.25:10; >>yi=interp1(x,y,xi); >> plot(x,y,’0’,xi,yi) 例 计算上例的函数插值 >>y2=interp1(x,y,x0); >>plot(x0,y2.’*m’)

  7. 7.1.3 Hermite插值 1.方法介绍 已知n个插值节点x1,x2,…，xn及其对应的函数值 y1,y2,…,yn和一阶导数值y‘1,y’2,…,y‘n。则计算插值区域内任意x的函数值y的Hermite插值公式：

  8. 2 MATLAB实现 hermite.m function y=hermite(x0,y0,y1,x) %hermite insert n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j- =i h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;

  9. a=1/(x0(i)-x0(j))+a; end end yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end 例 Hermite插值示例。对给定数据表如表7.2所示，试构造Hermite多项式，并给出sin0.34的近似值。 >> x0=[0.3 0.32 0.35]; >>y0=[0.29552 0.31457 0.34290];

  10. >>y1=[0.95534 0.94924 0.93937]; >>x=[0.3:0.005:0.35]; >>y=hermite(x0,y0,y1,x); >>plot(x,y) >>hermite (x0,y0,y1,0.34); >>y >>sin(0.34) >>y2=sin(x); >>hold on >>plot(x,y2,’--r’)

  11. 7.1.4 三次样条插值1 方法介绍设区间[a,b] 上给定有关划分 a=x0<x1<…<xn=b,S为[a,b] 上满足下列条件的函数. S∈C2 (a,b); S在每个子区间[xi,xi+1 ]上是三次多项式. 则称S关于划分的三次样条函数.常用的三次样条函数的边界条件有三种类型: I型 II型 其特殊情况为 此条件称为周期样条函数. Ⅲ型

  12. 2 MATLAB实现 Spline2.m function s=spline2(x0,y0,y21,y2n,x) n=length(x0); km=length(x); a(1)=-0.5; b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1))); for j=1:(n-1) h(j)=x0(j+1)-x0(j); end for j=2:(n-1) alpha(j)=h(j-1)/(h(j-1)+h(j));

  13. beta(j)=3*((1-alpha(j))*(y0(j)-y0(j-1))/h(j-1)+… alpha(j)*(y0(j+1)-y0(j))/h(j)); a(j)=-alpha(j)/(2+(1-alpha(j))*a(j-1)); b(j)=(beta(j)-(1-alpha(j))*b(j-1))/(2+(1-alpha(j))*a(j-1)); end m(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1)); for j=(n-1):-1:1 m(j)=a(j)*m(j+1)+b(j); end

  14. for k=1:km for j=1:(n-1) if((x(k)>=x0(j))&(x(k)<x0(j+1))) 1(k)=j; end end end for k=1:km sum=(3*(x0(1(k)+1)-x(k))^2/h(1(k))^2-… 2*(x0 (1(k)+1)-x(k))^3 /h(1(k))^3)*y0(1(k));

  15. sum=sum+(3*(x(k)-x0(1(k)))^2/h(1(k))^2-… 2*(x(k)-x0(1(k)))^3/h(1(k))^3/h(1(k))^3)*y0(1(k)+1); sum=sum+h(1(k))*((x0(1(k)+1)-x(k))^2/h(1(k))^2-…(x0(1(k) +1)-x(k))^3/h(1(k)^3)*m(1(k)); s(k)=sum-h(1(k))*((x(k)-x0(1(k)))^2/h(1(k))^2-… (x(k)-x0(1(k)))^3/h()1(k))^3)*m(1(k)+1); end

  16. 例 给定如下数据表所示的数据，试求三次样插值函数满足边界条件 s〞(28.7)= s〞(30)=0. 表 X 28.7 28 29 30 Y 4.1 4.3 4.1 3.0 >>x=[28.7 28 29 30]; >>y=[4.1 4.3 4.1 3.0]; >>x0=[28.7:0.15:30]; >>y=spline2(x,y,0,0,x0); >>plot(x0,y)

  17. 7.1.5 最小二乘法拟合 1 利用polyfit函数进行拟合 例：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) x1=[0.5:0.05:3.0]; y1=a(1)*x1.^2 +a(2)*x1 + a(3); plot(x,y,’*’,x1,y1,’-r’) 2 利用矩阵除法解决复杂型函数的拟合 例：用最小二乘法求一个形如y=a+bex的经 验公式 ，其中a、b为待定参数

  18. x=[0:0.4:2.0]; y=[1.2 1.3 1.45 1.6 2.0 2.5]; x1=exp(x) x1=[ones(6,1),x1’] %由y’=x1ab 得 ab=x1\y’ %其中 ab= x2=0:0.01:2.0; y2=ab(1)+ab(2)*exp(x2); plot(x,y,’*’,x2,y2,’r’) a b

  19. 7.1.6 快速Fourier变换简介 周期性的连续函数可以用三角多项式逼近。对 一组周期性的离散数据，这种逼近用快速Fourier 变换来实现。 设数据列为 x(n), n=1,2,…,N, 它对应的Fourier变换为 X(k), k=1,2,…,N 则有 X(k)= ∑ x(n)*exp(-j*2*pi*(k-1)*(n-1)/N) 1<= k <=N x(n)=(1/N) ∑X(k)*exp(j*2*pi*(k-1)*(n-1)/N) 1<= n <=N … N n=1 N n=1

  20. 命令函数及句法结构： fft(x) fft(x,n) fft(x,[],DIM)或fft(x,n,DIM) fft2 二维快速离散Fourier变换 调用形式(逆函数相同) fft2(x) fft2(x,MROWS,NCOLS) fftn n维快速离散Fourier变换 调用形式(逆函数相同) fftn(x) fftn(x,siz)

  21. ifft 快速离散Fourier逆变换 ifft2 二维快速离散Fourier逆变换 ifftn n维快速离散Fourier逆变换 例：给出一张记录 {xk }=[4 3 2 1 0 1 2 3],用FFT算法求{xk},离散频谱{Ck},其中k=0,1, …,7. 解: >>x=[4 3 2 1 0 1 2 3]; >>fft(x)

  22. 卷积运算演示 >> x=[4 3 2 1 0 1 2 3]; >> y=[1 1 1 1 1 1 1 1]; >>fx=fft(x,19) >>fy=fft(y,19) >>xy=conv(x,y) >>yx=ifft(fx.*fy) >>real(yx) 例：见demo/numerics/fast fourier transform

  23. 7.2 积分与微分 7.2.1 Newton-Cotes系列数值求值公式 若将积分区间[a,b]划分为n等份,步长 h=(b-a)/n,选取等距节点xk =a+kh,构造下面的插值型求公式. In =(b-a) 其中xk =a+kh,Ck(n)由Cotes系数表7.6给出. 当n=0时,即为矩形公式(稍做变化) 当n=1时,即为梯形公式 当n=2时,即为Simpson公式 当n=3时,即为Cotes公式

  24. 1矩形求积公式 y=cumsum(x) 对n维向量x， 返回n维向量y，y(i)是x的前i个分量之和。 例: >>x=[0 1 2;3 4 5]; >>cumsum(x,1) >>cumsum(x,2)

  25. 2 梯形求积公式 y=trapz(x) 对n维向量x，返回标量数值y，它的值是x的 各分量相隔为1单位时各小梯形面积之和。 例：求积分 e-0.5t sin(t+/6)dt 解：建立被积函数 function y=fun1(t) y=exp(-0.5*t).*sin(t+pi/6); ∫ 3π 0

  26. d=pi/1000; t=0:d:3*pi; nt=length(t); y=fun1(t); sc=cumsum(y)*d; scf=sc(nt) z=trapz(y)*d 3 自适应法－－Simpson法求积 Q=quad(‘F’,a,b) 求函数F(x) 从a到b的相对误差为1e-3的积分近似值 Q=quad(‘F’,a,b,tol) tol分别用于指定相对误差和绝对误差 例： ∫dx x 1 ————— x2+4 0

  27. function f=fun2(x) f=x./(4+x.^2) quad(‘fun2’,0,1) • 自适应法的Cotes求积公式 quad8 自适应法的Cote求积公式(高阶方法) Q=quad8(‘F’,a,b)使用Newton-Cotes法则求自适应递归法求函数F(X)从a到b的相对误差为1e-3的积分近似值,其中F为函数名. Q=quad8(‘F’,a,b,tol) tol分别用于指定相对误差和绝对误差

  28. 例：求积分 ∫e－x/2 dx 编制函数文件: fun.m function f=fun(x) f=exp(-x/2); >>quad8(‘fun’,1,3,1e-10)

  29. 7.2.2 Gauss求积公式 对任意的求积区间[a,b],通过变换x=(b- a)x/2+(a+b)/2,可转化到区间[-1，1]上。此时

  30. 例 求下面的积分值 先确定n的值。由高斯余项可得下式： 编制函数文件如下： gaussf.m function y=gaussf(x) y=cos(x); >>gauss1(0,1,6)

  31. 以n=3为例 gauss3.m function g=gauss3(fun,a,b,tol) if nargin<4 tol=1e-6;end c=(a+b)/2; x=[a,b,c]; y=feval(fun,x); lev=1; if any(imag(y)) g0=1e30; else

  32. g0=inf; end g=gaussstp1(fun,a,b,tol,lev,g0); function g=gaussstp1(fun,a,b,tol,lev,g0) levmax=10; if lev>levmax disp(‘Reearsion level limit reached in gauss’); g=g0; else c=(a+b)/2; g1=gaussstp2(fun,a,c); g2=gaussstp2(gun,c,b); g=g1+g2; if abs(g-g0)>tol*abs(g)

  33. g1=gaussstp2(fun,a,c,tol/2,lev+1,g1); g2=gaussstp2(fun,c,b,tol/2,lev+1,g2); g=g1+g2; end end function g=gaussstp2(fun,a,b) h=(b-a)/2; c=(a+b)/2; x=[h*(-0.7745967)+c,c,h*0.7745976+c]; f=feaval(fun,x); g=h*(0.555555556*(f(1)+f(3))+0.888888889*f(2)); >>gauss3(‘guassf’,0,1)

  34. 7.2.3 Romberg求积公式 1 Romberg加速法 rbgl.m function [s,n]=rbg1(a,b,eps) if nargin<3,eps=1e-6;end s=10; s0=0; k=2; t(1,1)=(b-a)*(f(a)+f(b))/2;

  35. while(abs(s-s0)>eps) h=(b-a)/2^(k-1); w=0; if(h~=0) for i=1:(2^(k-1)-1) w=w+f (a+i*h); end t(k,l)=h*(f(a)/2+w+f(b)/2); for l=2:k for i=1:(k-l+1) t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1); end end

  36. s=t(1,k); s0=(t(1,k-1)); k=k+1; n=k; else s=s0; n=-k; end end 例 用Romberg求积分公式计算积分

  37. 解：改编f.m函数如下： function f=f(x); f=x.^3; 计算积分，在MATLAB命令窗口中输入： rbg1(0,2) 2 改进Romberg的求积函数 rbg2.m

  38. function [s,eer]=rbg2(a,b,eps) if nargin<3, eps=1e-6;end m=1; t(1,1)=(b-a)*(f(a)+f(b))/2; r(1,1)=0; while((abs(t(1,m)-r(1,m))/2)>eps) c=0; for j=1:2^(m-1) c=c+f(a+(j-0.5)*(b-a)/2^(m-1)); end r(m,1)=(b-a)*c/2^(m-1); for j=2:m for k=1(m-j+1)

  39. r(k,j)=r(k+1,j-1)+(r(k+1,j-1)-r(k,j-1))/(4^(j-1)-1); end t(l,j)=r(l,j-1)+2*(4^(j-2)-1)*(t(1,j-1)-r(1,j-1))/(4^(j-1)-1); end end err=abs(t(1,m)-r(1,m))/2; s=t(1,m); 例 计算积分值

  40. 解：被积分的M文件f.m变为: function f=f(x) f=x.*sqrt(1+x.^2); 计算积分，在MATLAB命令窗口中输入： >>rbg1(0,3) 7.2.4 Mote-Carlo方法简介（积分多维） 求解半球体积 设半球体积为： f(x,y,z)= 其中

  41. 由于是半球，因此， 编制Mote-Carlo函数 mtc1.m function y=mtc1(fun,x,y,z,n) if nargin<5 n=1000;end n0=0; xh=x(2)-x(1); yh=y(2)-y(1); zh=z(2)-z(1); for i=1 a=xh*rand+x(1) b=yh*rand+y(1) c=zh*rand+z(1)

  42. f=feval(fun,a,b,c); if f<=0 n0=n0+1; end end y=n0/n*xh*yh*zh 球函数 f.m function f=f(x,y,z); f=x.^2+ y.^2+ z.^2-3; 积分计算：

  43. >>x=[-sqrt(3),sqrt(3)] >>y=[-sqrt(3),sqrt(3)] >>z=[0,sqrt(3)]; >>mtc1(‘f’,x,y,z,) 当增加随机点的个数时 mtc1(‘f’,x,y,z,1000) 7.2.6 微分和差分 1 数值微分和差分 diff(x) 对向量x ,其值为[ x(2)-x(1) x(3)-x(2) … x(n)-x(n-1)] diff(X) 对矩阵X，其值为矩阵列的差分[X(2：n,:)-X(1:n-1,:)]

  44. diff(X,n) 求n阶的差分值，如果n size(X,DIM), diff函数将先计算可能的连续差分值，直到size(X,DIM)=1. 然后diff 沿任意 n+1 维进行差分计算 diff(X,n,DIM) 求n阶的差分值，如果n size(X,DIM),函数返回空数组 2 符号微分和差分 diff(S) 对由findsym返回的自变量微分符号表达式S diff(S，’V’) 或diff(S，sym’V’) 对自变量v微分符号表达式S diff(S，n) 对正整数n，函数微分S表达式n次 diff(S,n,’V’) 和diff(S,’V’ , n）两种形式都可以被识别

  45. 例 x=sym(‘x’); diff(sin(x^2)) 3 梯度函数 gradient 近似梯度函数 [fx,fy]=gradient(F) 返回矩阵F的数值梯度，fx相当于 dF/dx, 为x 方向的差分值。fy相当于dF/dx,为y方向的差分值。各个方向的点间隔设为1。当F为向量分值，DF=gradient(F) 为一维梯度

  46. 7.3 求一元高次方程与线性方程组 an xn+an-1 xn-1+…+a0=0 p=[an ,an-1 , …,a0 ] x0=roots(p) n元线性方程组可表达为 AX＝b ％A为nxn阶系数矩阵 ％b为nx1阶常数矩阵 ％ X为nx1阶变量矩阵 X＝A\b

  47. 例 求解下列方程组 解：在MATLAB命令窗中输入： a=[0.4096 0.1234 0.3678 0.2943 0.2246 0.3872 0.4015 0.1129 0.3645 0.1920 0.3781 0.0643 0.1784 0.4002 0.2786 0.3927]； b=[0.4043 0.1550 0.4240 –0.2557]’； x=a\b

  48. 7.3.2 迭代解法的几种形式 1 Jacobi迭代法 方程组Ax=b,其中A∈ ，b ∈ , 且A为非奇异，则A可写成A=D-L-U。其中，D= diag , 而-L、-U分别为A的严格下三角部分（不包括对角线元素） 则

  49. 2 Gauss-Seidel(G-S)迭代法 3. SOR迭代法 4.两步迭代法 当线性方程系数矩阵为对称正定时，可用一种特殊的迭代法来解决，其迭代公式如下：

  50. MATLAB的实现 jacobi.m D=diag(diag(a)); U=-triu(a,-1); L=-triu(a,-1); B=D\(L+U); f=D\b; y=B*x0+f;n=1; while norm(y-x0)>1.0e-6 x0=y; y =B*x0+f;n=n+1; end y n

More Related