《计算方法的matlab实现》由会员分享，可在线阅读，更多相关《计算方法的matlab实现（23页珍藏版）》请在金锄头文库上搜索。
1、误差基础% ex01.mclear;digit = 3;S0 = log(6)-log(5); S0 = vpa(S0,digit);S(1) = 1 - 5*S0;S0 = vpa(S(1),digit);for n = 2 : 8S(n) = 1/n - 5*S(n-1);S(n) = vpa(S(n),digit);endS=double(S);fprintf( S0 = %11fn,double(S0);for n=1:8fprintf( S%1d = %11fn,n,S(n);end_% ex02.mclear;digit = 3;S(8) = (1/45 + 1/54)/2; S(
2、8) = vpa(S(8),digit);for n = 7 : -1: 1 S(n) = 1/(5*(n+1) - S(n+1)/5;S(n) = vpa(S(n),digit);endS0 = 1/5 - S(1)/5;S0 = vpa(S0,digit);S=double(S);fprintf( S0 = %11fn,double(S0);for n=1:8fprintf( S%1d = %11fn,n,S(n);end_%ex03.mclear;disp( 用单精度计算 108+1);result = single(1e8 + 1);fprintf( 108+1 = %fn, resu
3、lt);fprintf(n press any key to continue . n); pause;disp( 计算 1 + 2 + . + 40 + 109)result = 0;for k = 1 : 40result = single(result + k);endresult = single(result+1e9);fprintf( 1 + 2 + . + 40 + 109 = %fn, result);fprintf(n press any key to continue . n); pause;disp( 计算 109 + 1 + 2 + . + 40)result = 1e
4、9;for k = 1 : 40result = single(result + k);endfprintf( 108 + 1 + 2 + . + 40 = %fn, result);_lagrange 插值% textbook page 35function yh = lagrange(x,y,xh)n = length(x);m = length(xh);x = reshape(x,n,1); % x = x(:);y = reshape(y,n,1); % y = y(:);xh = reshape(xh,m,1); % xh = xh(:);yh = zeros(m,1);c1 = o
5、nes(1,n-1);c2 = ones(m,1);for i=1:n,xp = x(1:i-1,i+1:n);yh = yh + y(i) * prod(xh*c1-c2*xp)./(c2*(x(i)*c1-xp),2);end_% Lagrange interpolation 1function yh = mylagrange1(x,y,xh)n = length(x)-1;x = reshape(x,1,n+1); y = reshape(y,1,n+1); yh = 0;for j = 0 : n xp = x(1:j, j+2:n+1);yh = yh + y(j+1) * prod
6、(xh - xp)./(x(j+1) - xp);end_% Lagrange interpolation 2function yh = mylagrange1(x,y,xh)% xh can be a vectorn = length(x)-1;m = length(xh);x = reshape(x,1,n+1);y = reshape(y,1,n+1);xh = reshape(xh,m,1);yh = zeros(m,1);for j = 0 : n xp = x(1:j, j+2:n+1);tmp = (xh*ones(1,n)-ones(m,1)*xp)./(ones(m,1)*(
7、x(j+1) - xp);yh = yh + y(j+1)*prod(tmp,2);end_newton 插值function p,q = chashang(x,y)n = length(x);x = reshape(x,n,1);p(:,1) = x;p(:,2) = reshape(y,n,1);for j = 3: n+1p(1:n+2-j,j) = diff(p(1:n+3-j,j-1) ./ (x(j-1:n)-x(1:n+2-j);endq = p(1,2:n+1);_function yh = newton(x,y,xh) p,q = chashang(x,y); n = len
8、gth(x); m = length(xh); x = reshape(x,1,n);xh = reshape(xh,m,1);xh = xh*ones(1,n) - ones(m,1)*x; xh = xh; yh = y(1)*ones(1,m); for i = 2 : n yh = yh + q(i)*prod(xh(1:i-1,:),1); end _-三次样条插值和 Matlab 插值函数介绍% textbook page 47clear; close all; clc;xh = -5:1/4:5;yt = 1./(1+xh.2); % true solutionplot(xh,y
9、t,k,linewidth,2);axis(-5 5, -0.5 2.5);hold on;color=r,b,g;nn = 2 5 10;for k = 1:length(nn)n = nn(k);x = -5:10/n:5;y = 1./(1+x.2);yh = lagrange(x,y,xh); % Lagrange interpolationpauseplot(xh,yh,colork,linewidth,2);endhold off;_% 比较 f(x)=1/(1+x2) 的各种插值函数clear; close all; a = -5; b = 5; % 给定插值区间n = 10;
10、h = (b-a)/n; % 计算步长（等距节点个数为 n+1）x = a:(b-a)/n:b; y = 1./(1+x.2); % 计算节点和函数值xh = a:(b-a)/40:b; % 插值点% 首先画出原函数图形yh1 = 1./(1+xh.2); % true solutionplot(xh,yh1,k.-,linewidth,1.5);axis(a, b, -0.5 2.5);hold on;% 10 次 Lagrnage 插值yh2 = lagrange(x,y,xh); % 由于 Lagrange 插值不稳定，系统不提供该函数plot(xh,yh2,b+-,linewidth,
11、1.5);% 3 次样条插值df = diff(1/(1+x2),x);df0 = subs(df,x,-5); dfn = subs(df,x,5); % 计算第一类边界条件yh3 = spline(x,df0 y dfn, xh);plot(xh, yh3, ro-, linewidth,1.5);% 线性插值yh4 = interp1(x,y,xh);plot(xh, yh4, gd-, linewidth,1.5);legend(原函数图形,10 次 Lagrange插值, .三次样条插值,线性插值);plot(x,y,s,linewidth,2); % 画出点 (xi,f(xi)ho
12、ld off;_% myspline02.mclear;a = -5; b = 5; n = 40; h = (b-a)/n;xh = a:h:b;syms x fx = (x2+1)(-1);dfx = diff(fx);dfa = subs(dfx,x,a);dfb = subs(dfx,x,b);yh = subs(fx,x,xh);mu = ones(1,n-1)/2;lambda = ones(1,n-1)/2;p,q = chashang(xh,yh);d = zeros(n+1,1);d(2:n) = 6*p(1:n-1,4);d(1) = 6*(p(1,3) - dfa)/h;
13、d(n+1) = 6*(dfb - p(n,3)/h;A = diag(2*ones(1,n+1) + diag(mu,1,-1) + diag(1,lambda,1);M = Ad;a3 = (M(2:n+1)-M(1:n)/(6*h);a2 = M(1:n)/2;a1 = (yh(2:n+1) - yh(1:n)/h - h*(M(2:n+1) +2*M(1:n)/6;a0 = yh(1:n);pp = spline(xh,dfa yh dfb);breaks,coefs,npolys,ncoefs,dim = unmkpp(pp);fprintf(n err_a3=%f, err_a2=
14、%f, err_a1=%f, err_a0=%fn, . norm(a3-coefs(:,1), norm(a2-coefs(:,2), .norm(a1-coefs(:,3), norm(a0-coefs(:,4);% plot the figuresnp = 3;n0 = np*n; h = (b-a)/n0;xh0 = a:h:b; len = length(xh0);xh =xh(1:n);for k = 1 : npindex = k:np:len;index = index(1:n);xht = xh0(index); yht = a3.*(xht-xh).3 + a2.*(xht-xh).2 + a1.*(xht-xh) + a0;yh0(index) = yht;endyh0(n0+1) = yh(n+1);yh1 = subs(fx,x,xh0);plot(xh0,yh1,b.-,xh0,yh0,ro-)_
