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

Matlab在计算方法中的应用1.插值与拟合2.积分与微分3.求解线型方程组4.求解非线性方程组5.常微分方程的解法1.插值与拟合1.Lagrange插值%lagrangeinsertfunctiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endendnknkjjjkjkxxxxyxy11)()(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的函数值y2.分段线型插值：所谓分段线型插值就是通过插值点用折线段连接起来逼近原曲线1.yi=interp1(x,y,xi)对节点向量（x,y）插值，求xi对应的yi值2.yi=interp1(y,xi)默认x=1:n，n为向量y的长度值3.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、interpn012345678910-1-0.8-0.6-0.4-0.200.20.40.60.813.Hermite插值要求插值点上函数值和导数值都相等xi,yi,yi’分别为插值节点、对应函数值和对应一阶倒数值。自编程序函数：y=hermite(x0,y0,y1,0.34);nijjnijijiijijiniiiiiiixxaxxxxhyyyaxxhxy1121'1;)(])2)([()(其中：functiony=hermite(x0,y0,y1,x)%hermiteinsertn=length(x0);m=length(x);fork=1:myy=0.0;fori=1:nh=1.0;a=0.0;forj=1:nifj~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end4.三次样条插值设区间[a,b]上给定的有关划分a=x0x1…xn=b,S为[a,b]区间上满足下面条件的函数：•S在[a,b]上二阶导数连续•S在每个插值子区间[xi,xi+1]上是三阶多项式则称S为关于划分的有关三次样条函数。常用的三次样条函数的边界条件有三种类型•I型，S’(x0)=f’0,S’(xn)=f’n;•II型，S’’(x0)=f’’0,S’’(xn)=f’’n;特殊情况为都等于0；•III型，Sj(x0)=Sj(xn)，j=0,1,2,…;周期样条函数。自编II型程序函数如下：s=spline2(x0,y0,y21,y2n,x)functions=spline2(x0,y0,y21,y2n,x)%s=spline2(x0,y0,y21,y2n,x)%x0,y0areexistedpoints,xareinsertpoints,%y21,y2nareseconddirivitivenumbersgivenn=length(x0);km=length(x);a(1)=-0.5;b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1)));forj=1:(n-1)h(j)=x0(j+1)-x0(j);endforj=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));endm(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1));forj=(n-1):-1;1m(j)=a(j)*m(j+1)+b(j);endfork=1:kmforj=1:(n-1)if((x(k)=x0(j))&(x(k)x0(j+1)))l(k)=j;endendendfork=1:kmsum=(3*(x0(l(k)+1)-x(k))^2/h(l(k))^2-...2*(x0(l(k)+1)-x(k))^3/h(l(k))^3)*y0(l(k));sum=sum+(3*(x(k)-x0(l(k)))^2/h(l(k))^2-...2*(x(k)-x0(l(k)))^3/h(l(k))^3)*y0(l(k)+1);sum=sum+h(l(k))*((x0(l(k)+1)-x(k))^2/h(l(k))^2-...(x0(l(k)+1)-x(k))^3/h(l(k))^3)*m(l(k));s(k)=sum-h(l(k))*((x(k)-x0(l(k)))^2/h(l(k))^2-...(x(k)-x0(l(k)))^3/h(l(k))^3)*m(l(k)+1);end4.最小二乘拟合1.利用polyfit进行多项式拟合x=[0.51.01.52.02.53.0];y=[1.752.453.814.807.008.60];a=polyfit(x,y,2)a=0.56140.82871.1560x1=[0.5:0.05:3.0];y1=a(3)+a(2)*x1+a(1)*x1.*x1;plot(x1,y1,'-r')holdonplot(x,y,'*')0.511.522.531234567892.利用常用的矩阵除法解决复杂型函数的拟合例：用最小二乘法求形如y=a+bx2的经验公式，使它与下表数据拟合xi1925313844yi19.032.349.073.397.8x=[1925313844];y=[1932.34973.397.8];x1=x.^2;x1=[ones(5,1),x1']x1=1361162519611144411936ab=x1\y'ab=0.97260.0500x0=[19:0.2:44];y0=ab(1)+ab(2)*x0.^2;clfplot(x,y,'o')holdonplot(x0,y0,'-r')152025303540451020304050607080901002.积分与微分1.Newton-Cotes系列数值求积公式1.矩形求积公式cumsum(X)2.梯形求积公式trapz(X,Y)3.自适应simpson法求积quad(‘F’,a,b,…)4.自适应的cotes法求积公式quad8(‘F’,a,b,…)2.Gauss求积公式3.Romberg求积公式4.Monte－Carlo方法以上均可自编程序完成。5.符号积分intsymsum6.微分和差分1.数值微分与差分diff（X,N,DIM）2.符号微分与差分diff（S，’v’，n）3.梯度函数[fx,fy]=gradient(F,HX,HY)4.多元函数的导数jacobian(f,v)3.求解线型方程组一般分为两种•直接法：通过矩阵的变形、消去直接求解，主要用于低阶稠密矩阵•叠代法：利用某种极限过程去逐渐逼近方程组精确解，主要用于大型稀疏矩阵1.直接法：1.矩阵除法：x=a\b2.线性方程组直接求解分析1.LU分解：[l,u]=lu(a)2.Cholesky分解：l=chol(a)3.奇异值分解：[U,S,V]=SVD(X)4.上三角变换：triu5.对角变换：diag6.下三角变换：tril2.跌代解法的几种形式1.Jacobi跌代法2.gauss-seidel跌代法3.SOR（逐次超松弛跌代法）4.两步跌代法均可自己编程完成3.线性方程组的解析解法linsolvesolvevpa4.求解非线性方程组1.非线性方程的解法1.二分法2.不动点叠代法3.Newton叠代（切线叠代）法4.割线法可自行编制函数2.方程组解法1.不动点跌代2.Newton法3.broyden法（秩1的拟newton法）3.非线性方程（组）的解析解法fsolve(‘fc’,x0)5.常微分方程的解法1.欧拉方法1.简易欧拉法2.改进欧拉法需自编程序完成•Runge-Kutta方法ODE解函数：ode23,ode45,ode113,ode15s,ode23s参数选择函数：odeset,odeget输出函数：odeplot,odephas2,odephas3,odeprintODE范例：orbt2ode,rigidobe,vdpode1)二三阶R-K函数（低阶方法）[T,Y]=ode23(‘F’,TSPAN,Y0,…)F为求解微分方程，TSPAN为微分方程积分限，Y0初始条件2)四五阶P-K函数（中阶方法）[T,Y]=ode45(‘F’,TSPAN,Y0,…)F为求解微分方程，TSPAN为微分方程积分限，Y0初始条件3.常微分方程的解析解dsolve前面已作介绍