Matlab在计算方法中的应用.ppt

Matlab在计算方法中的应用,插值与拟合 积分与微分 求解线型方程组 求解非线性方程组 常微分方程的解法,1.插值与拟合,Lagrange插值 %lagrange insert 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 if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end end, x=0.4:0.1:0.8; y=-0.916 -0.693 -0.510 -0.356 -0.223; lagrange(x,y,0.54) ans = -0.6161,对给定n个插值节点x1,x2,xn及对应函数值y1,y2,yn，利用(n-1)次lagrange插值多项式公式，求得插值区间任意x的函数值y,分段线型插值： 所谓分段线型插值就是通过插值点用折线段连接起来逼近原曲线 yi=interp1(x,y,xi) 对节点向量（x,y）插值，求xi对应的yi值 yi=interp1(y,xi) 默认x=1:n，n为向量y的长度值 yi=interp1(x,y,xi，method) method指定插值的算法，默认为线型算法，可取值为：nearest线性最近项插值；linear线性插值；spline立方样条插值；cubic立方插值。, x=0:0.1:10; y=sin(x); xi=0:0.25:10; yi=interp1(x,y,xi); plot(x,y,o,xi,yi),同类函数：interp1q、interpft、spline、interp2、interp3、interpn,Hermite插值 要求插值点上函数值和导数值都相等 xi, yi, yi分别为插值节点、对应函数值和对应一阶倒数值。 自编程序函数：y=hermite(x0,y0,y1,0.34);,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; 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,三次样条插值 设区间a,b上给定的有关划分a=x0x1xn=b, S为a,b区间上满足下面条件的函数： S在a,b 上二阶导数连续 S在每个插值子区间xi,xi+1上是三阶多项式 则称S为关于划分的有关三次样条函数。 常用的三次样条函数的边界条件有三种类型 I型，S(x0)=f0, S(xn)=fn; II型，S(x0)=f0, S(xn)=fn;特殊情况为都等于0； III型， Sj(x0)= Sj(xn)，j=0,1,2,; 周期样条函数。 自编II型程序函数如下：s=spline2(x0,y0,y21,y2n,x),function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points, x are insert points, %y21,y2n are second dirivitive numbers given 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); 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,for k=1:km for j=1:(n-1) if(x(k)=x0(j) end,最小二乘拟合 利用polyfit进行多项式拟合 x=0.5 1.0 1.5 2.0 2.5 3.0; y=1.75 2.45 3.81 4.80 7.00 8.60; a=polyfit(x,y,2) a = 0.5614 0.8287 1.1560 x1=0.5:0.05:3.0; y1=a(3)+a(2)*x1+a(1)*x1.*x1; plot(x1,y1,-r) hold on plot(x,y,*),利用常用的矩阵除法解决复杂型函数的拟合 例：用最小二乘法求形如y=a+bx2的经验公式，使它与下表数据拟合, x=19 25 31 38 44; y=19 32.3 49 73.3 97.8; x1=x.2; x1=ones(5,1),x1 x1 = 1 361 1 625 1 961 1 1444 1 1936, ab=x1y ab = 0.9726 0.0500 x0=19:0.2:44; y0=ab(1)+ab(2)*x0.2; clf plot(x,y,o) hold on plot(x0,y0,-r),2.积分与微分,Newton-Cotes系列数值求积公式 矩形求积公式 cumsum(X) 梯形求积公式 trapz(X,Y) 自适应simpson法求积 quad(F,a,b,) 自适应的cotes法求积公式 quad8(F,a,b,),Gauss求积公式 Romberg求积公式 MonteCarlo方法 以上均可自编程序完成。 符号积分 int symsum,微分和差分 数值微分与差分 diff（X,N,DIM） 符号微分与差分 diff（S，v，n） 梯度函数 fx,fy=gradient(F,HX,HY) 多元函数的导数 jacobian(f,v),3.求解线型方程组,一般分为两种 直接法：通过矩阵的变形、消去直接求解，主要用于低阶稠密矩阵 叠代法：利用某种极限过程去逐渐逼近方程组精确解，主要用于大型稀疏矩阵,直接法： 矩阵除法：x=ab 线性方程组直接求解分析 LU分解： l,u=lu(a) Cholesky分解：l=chol(a) 奇异值分解：U,S,V = SVD(X) 上三角变换：triu 对角变换：diag 下三角变换：tril,跌代解法的几种形式 Jacobi跌代法 gauss-seidel跌代法 SOR（逐次超松弛跌代法） 两步跌代法 均可自己编程完成 线性方程组的解析解法 linsolve solve vpa,4.求解非线性方程组,非线性方程的解法 二分法 不动点叠代法 Newton叠代（切线叠代）法 割线法 可自行编制函数,方程组解法 不动点跌代 Newton法 broyden法（秩1的拟newton法） 非线性方程（组）的解析解法 fsolve(fc,x0),5.常微分方程的解法,欧拉方法 简易欧拉法 改进欧拉法 需自编程序完成 Runge-Kutta方法 ODE解函数：ode23, ode45, ode 113, ode15s, ode23s 参数选择函数：odeset, odeget 输出函数：odeplot, odephas2, odephas3, odeprint ODE范例：orbt2ode,