数值分析教学课件：Ch10-数值算法实现

Ch10.数值算法实现数值算法实现1.线性方程组解法线性方程组解法 1、三角形线性方程组解法、三角形线性方程组解法 以上三角形线性方程组 为例。回代：iinijjijiinnnnaxabxabx1,1,2,2,1nniAxb%文件 uptri.mfunction u=uptri(a,b)n=size(a,1);x(n)=b(n)/a(n,n);for i=n-1:-1:1 s=0;for j=i+1:n s=s+a(i,j)*x(j);end x(i)=(b(i)s)/a(i,i);endu=x;求解求解Axbnijjijxas12、顺序Gauss消去法 （1）消去过程：第 k 步，计算（2）回代过程：()()(1)()()(1)()(),1,1,;1,1,kkikikkkkkkijijikkjkkkiiikkaaiknaaaikn jknbbbikn)1,2,1(nk)(1)()()()(,iiinijjiijiiinnnnnnaxabxabx1,2,2,1nni%文件 gauss.mfunction u=gauss(a,b)n=length(b);for k=1:n 1 for i=k+1:n mult=a(i,k)/a(k,k);for j=k+1:n%if abs(a(k,k)1e6 a(i,j)=a(i,j)mult*a(k,j);%else%disp(顺序Gauss消去法失败);%pause；%exit；%end end可去掉%，判断主元是否为0 b(i)=b(i)mult*b(k);endendx(n)=b(n)/a(n,n);for i=n-1:-1:1 s=0;for j=i+1:n s=s+a(i,j)*x(j);end x(i)=(b(i)s)/a(i,i);endu=x;回代回代例：%主文件main.ma=6,-2,2,4;12,-8,6,10;3,-13,9,3;-6,4,1,-18;b=16,26,-19,-34;x=gauss(a,b);disp(方程组解为方程组解为:);x则有:main方程组解为方程组解为:x=3 1 -2 13、Jacobi迭代法迭代法 bAx ADL U000000022311312213231212211nnnnnnaaaaaaaaaaaaa,)(1)(1)1(bDxULDxkk0,1,2,k%jacobi.mfunction y=jacobi(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(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;endyn例:P249 例7.21 a=10,-1,0;-1,10,-2;0,-2,10;b=9;7;6;jacobi(a,b,0;0;0)y=0.9958 0.9579 0.7916n=11解向量解向量迭代次数迭代次数2.方程求根方程求根1、二分法、二分法%erfen.mfunction y=erfen(fun,a,b,esp)if nargin 4 esp=1e 4;endif feval(fun,a)*feval(fun,b)esp if feval(fun,a)*feval(fun,c)0 b=c;c=(a+b)/2;elseif feval(fun,c)*feval(fun,b)erfen(fc,0,10,0.05)n=56ans=1.6180先编写函数先编写函数2、不动点迭代法迭代法%iterate.mfunction y=iterate(x)x1=g(x);n=1;while (abs(x1 x)=1.0e 6)&(n iterate(3)x1=3.7331n=22初值不同初值不同,迭代步数会不同迭代步数会不同迭代函数迭代函数3、牛顿迭代法迭代法%newton.mfunction y=newton(x0)x1=x0 fc(x0)/df(x0);n=1;while(abs(x1 x0)=1.0e 6)&(n newton(0)newton(10)x1=x1=-0.4590 3.7331n=n=6 12比不动点比不动点迭代快迭代快 习题三习题三1、编写下三角形方程组前推法的函数；2、编写Jacobi迭代分量形式的函数；3、利用牛顿迭代法 求解计算方法P114 习题4.7，写出各函数和命令文件。下下 次次 见见3.插值方法插值方法 以牛顿插值为例。以牛顿插值为例。已知函数 在互异点 上的值为 。n n次次Newton插值多项式插值多项式 为：nnnxxxfxxxxxxxxxfxxxxxxfxxxfxN,)()(,)(,)()()(10110210101000)(xf,10 xx),(),(10 xfxf定义定义 设函数 在互异点 上的值为 。定义：1）一阶均差一阶均差为2）二阶均差二阶均差为3）递推下去，k 阶均差阶均差为,10 xx),(),(10 xfxfijijxxxfxfjixxf)()(,ikjikjxxxxfxxfkjixxxf,011021,10,xxxxxfxxxfkkkkxxxf)(xf例：已知 f(x)=sqrt(2x)+ln(x)节点 x1=0.5,x2=1,x3=1.5,x4=2.0 求 Newton 插值多项式，并计算其在 x=0.75,x=1.75处的值。比较与原函数的误差。%f.mfunction y=f(x)y=sqrt(2*x)+log(x);%j0.mfunction r=j0(x)r=f(x);%j1.mfunction r=j1(x1,x2)r=(j0(x1)-j0(x2)/(x1-x2);原函数原函数0 阶均差阶均差(原函数原函数)1 阶均差阶均差%j2.mfunction r=j2(x1,x2,x3)r=(j1(x1,x2)-j1(x2,x3)/(x1-x3);%j3.mfunction r=j3(x1,x2,x3,x4)r=(j2(x1,x2,x3)-j2(x2,x3,x4)/(x1-x4);2 阶均差阶均差3 阶均差阶均差%newton.mfunction r=newton(t)x=0.5,1.0,1.5,2.0 ;r=j0(x(1)+j1(x(1),x(2).*(t-x(1)+j2(x(1),x(2),x(3).*(t-x(1).*(t-x(2)+j3(x(1),x(2),x(3),x(4).*(t-x(1).*(t-x(2).*(t-x(3);插值函数插值函数)()(,)(,)(,3214321213211211xxxxxxxxxxfxxxxxxxfxxxxfxfy%main.mdisp(x=0.75处值与误差);t=0.75;l=newton(t)e=abs(f(t)-l)disp(x=1.75处值与误差);t=1.75;l=newton(t)e=abs(f(t)-l)row=0.5:0.1:2.0;y=f(row);N=newton(row);plot(row,y,b-,row,N,r-);legend(原函数,插值函数);mainx=0.75处值与误差l=0.9221e=0.0150 x=1.75处值与误差l=2.4228e=0.00774.曲线拟合的最小二乘法曲线拟合的最小二乘法 1、多项式拟合多项式拟合2、一般情形、一般情形1、多项式拟多项式拟合合 已知(xi,yi)(i=1,m)polyfit(x,y,n)x=x1,x2,xm y=y1,y2,ym 多项式次数多项式次数例:用三次多项式在 0,5 上拟合 ex x=0:0.1:5;y=exp(x);p=polyfit(x,y,3);s=polyval(p,x);plot(x,y,b*,x,s,“r-)legend(原曲线,拟合曲线)axis(0,5,0,52 )拟合函数2、一般情形、一般情形对应于解超定方程组 用矩阵除法实现例：用 y=a+bx2 拟合给定的一组数据：),2,1(),(niyxii解：使解：使 最小。最小。对应的对应的超定方程组为：超定方程组为：即即Matlab实现：实现：u=A y (见见 P228 例例7.8)niiiniiybxa12212)(),.,2,1(2niybxaiinnyyybaxxx2022221111记为记为 Au=y5.数值积分数值积分 复化梯形公式的实现:trapz(y)*h 即y=f(x0),f(xn)nnkkkThxfxf*)()(21101例:P233 例7.10%fun.mfunction y=fun(t)y=exp(-0.5*t).*sin(t+pi/6);d=pi/1000;t=0:d:3*pi;y=fun(t);z=trapz(y)*dz=0.90084027660688被积函数被积函数5.常微分方程数值解法常微分方程数值解法 以以EulerEuler法为例。法为例。例例:取步长取步长 h=0.2h=0.21)0()50(1yxxyy%f.mfunction r=f(x,y)r=y-x+1;%orig.mfunction r=orig(x,y)r=x+exp(x);右端函数右端函数准确解准确解%euler.m Euler方法的函数function X,Y =euler(a,b,h,c)%c=y0n=round(b-a)/h)+1;%计算点的个数x=zeros(n,1);%设定 x 维数y=zeros(n,1);%设定 y 维数x(1)=a;y(1)=c;%初始条件str=sprintf(x0=%g,y0=%g n,x(1),y(1);disp(str);%显示初始条件%对每个点用Euler公式计算for i=1:n-1 y(i+1)=y(i)+h*f(x(i),y(i);x(i+1)=x(i)+h;str=sprintf(%d x=%g y=%g e=%g,i,x(i+1),y(i+1),abs(y(i+1)-orig(x(i+1);disp(str);%显示计算结果与误差endX=x;Y=y;%返回计算结果%主文件 main.mclear;%清除内存变量a=0;b=5;h=0.2;c=1;%输入参数值 X,Y =euler(a,b,h,c);%调用Euler函数%绘图plot(X,orig(X),b-,X,Y,r-);legend(精确解,数值解);mainx0=0,y0=1 1 x=0.2 y=1.4 e=0.02140282 x=0.4 y=1.84 e=0.051824725 x=5 y=100.396 e=53.0169 谢谢 谢谢