Advanced Mathematics and Mechanics Applications Using MATLAB phần 3 pptx

2: %

3: % [area,leng,X,Y,closed]=curvprop(x,y,doplot)

4: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

5: % This function passes a cubic spline curve through

6: % a set of data values and computes the enclosed

7: % area, the curve length, and a set of points on

8: % the curve.

9: % 10: % x,y - data vectors defining the curve.

11: % doplot - plot the curve if a value is input for

12: % doplot. Otherwise, no plot is made. 13: % area - the enclosed area is computed. This

14: % parameter is valid only if the curve

15: % is closed and the boundary is traversed 16: % in counterclockwise. For a curve, the

17: % area agrees with a curve closed using 18: % a line from the last point to the

19: % origin, and a line from the origin to

20: % the first point. 21: % leng - length of the curve

22: % X,Y - set of points on the curve. The output

23: % intersperses three points on each segment 24: % between the starting data values.

25: % closed - equals one for a closed curve. Equals zero 26: % for an open curve.

27: %

28:

29: % For default test data, choose an ellipse with

30: % semi-diameters of 2 and 1.

31: if nargin==0 32: m=21; th=linspace(0,2*pi,m);

33: x=2*cos(th); y=sin(th); x(m)=2; y(m)=0; 34: end

35:

36: % Use complex data coordinates 37: z=x(:)+i*y(:); n=length(z); t=(1:n)’;

38: chord=sum(abs(diff(z))); d=abs(z(n)-z(1));

39:

40: % Use a periodic spline if the curve is closed

41: if d < (chord/1e6) 42: closed=1; z(n)=z(1); endc=5;

43: zp=spterp(t,z,1,t,endc);

44:

45: % Use the not-a-knot end condition for open curve

46: else

© 2003 by CRC Press LLC

47: closed=0; endc=1; zp=spterp(t,z,1,t,endc);

48: end

49:

50: % Compute length and area

51: % plot(abs(zp)),shg,pause 52: leng=spterp(t,abs(zp),3,n,1);

53: area=spterp(t,1/2*imag(conj(z).*zp),3,n,1);

54: Z=spterp(t,z,0,1:1/4:n,endc); 55: X=real(Z); Y=imag(Z);

56: if nargin>2

57: plot(X,Y,’-’,x,y,’.’), axis equal 58: xlabel(’x axis’), ylabel(’y axis’)

59: title(’SPLINE CURVE’), shg

60: end 61:

62: %============================================ 63:

64: function [v,c]=spterp(xd,yd,id,x,endv,c)

65: % 66: % [v,c]=spterp(xd,yd,id,x,endv,c)

67: % ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

68: % 69: % This function performs cubic spline interpo￾70: % lation. Values of y(x),y’(x),y’’(x) or the 71: % integral(y(t)*dt, xd(1)..x) are obtained.

72: % Five types of end conditions are provided.

73: % 74: % xd, yd - data vectors with xd arranged in

75: % ascending order.

76: % id - id equals 0,1,2,3 to compute y(x), 77: % y’(x), integral(y(t)*dt,t=xd(1)..x),

78: % respectively. 79: % v - values of the function, first deriva￾80: % tive, second derivative, or integral

81: % from xd(1) to x 82: % c - the coefficients defining the spline

83: % curve. If these values are input from

84: % an earlier computation, then they 85: % are not recomputed.

86: % endv - vector giving the end conditions in 87: % one of the following five forms:

88: % endv=1 or endv omitted makes

89: % c(2) and c(n-1) zero 90: % endv=[2,left_end_slope,...

91: % right_end_slope] to impose slope

© 2003 by CRC Press LLC

92: % values at each end

93: % endv=[3,left_end_slope] imposes the

94: % left end slope value and makes 95: % c(n-1) zero

96: % endv=[4,right_end_slope] imposes the 97: % right end slope value and makes

98: % c(2) zero

99: % endv=5 defines a periodic spline by 100: % making y,y’,y" match at both ends

101:

102: if nargin<5 | isempty(endv), endv=1; end 103: n=length(xd); sx=size(x); x=x(:); X=x-xd(1);

104:

105: if nargin<6, c=spcof(xd,yd,endv); end 106:

107: C=c(1:n); s1=c(n+1); m1=c(n+2); X=x-xd(1); 108:

109: if id==0 % y(x)

110: v=yd(1)+s1*X+m1/2*X.*X+... 111: powermat(x,xd,3)*C/6;

112: elseif id==1 % y’(x)

113: v=s1+m1*X+powermat(x,xd,2)*C/2; 114: elseif id==2 % y’’(x)

115: v=m1+powermat(x,xd,1)*C; 116: else % integral(y(t)*dt, t=xd(1)..x)

117: v=yd(1)*X+s1/2*X.*X+m1/6*X.^3+...

118: powermat(x,xd,4)*C/24; 119: end

120: v=reshape(v,sx);

121:

122: %==============================================

123:

124: function c=spcof(x,y,endv)

125: %

126: % c=spcof(x,y,endv) 127: % ~~~~~~~~~~~~~~~~

128: % This function determines spline interpolation

129: % coefficients consisting of the support 130: % reactions concatenated with y’ and y’’ at

131: % the left end. 132: % x,y - data vectors of interplation points.

133: % Denote n as the length of x.

134: % endv - vector of data for end conditions 135: % described in function spterp.

136: %

© 2003 by CRC Press LLC

137: % c - a vector [c(1);...;c(n+2)] where the

138: % first n components are support

139: % reactions and the last two are 140: % values of y’(x(1)) and y’’(x(1)).

141:

142: if nargin<3, endv=1; end

143: x=x(:); y=y(:); n=length(x); u=x(2:n)-x(1);

144: a=zeros(n+2,n+2); a(1,1:n)=1; 145: a(2:n,:)=[powermat(x(2:n),x,3)/6,u,u.*u/2];

146: b=zeros(n+2,1); b(2:n)=y(2:n)-y(1);

147: if endv(1)==1 % Force, force condition 148: a(n+1,2)=1; a(n+2,n-1)=1;

149: elseif endv(1)==2 % Slope, slope condition

150: b(n+1)=endv(2); a(n+1,n+1)=1; 151: b(n+2)=endv(3); a(n+2,:)=...

152: [((x(n)-x’).^2)/2,1,x(n)-x(1)]; 153: elseif endv(1)==3 % Slope, force condition

154: b(n+1)=endv(2); a(n+1,n+1)=1; a(n+2,n-1)=1;

155: elseif endv(1)==4 % Force, slope condition 156: a(n+1,2)=1; b(n+2)=endv(2);

157: a(n+2,:)=[((x(n)-x’).^2)/2,1,x(n)-x(1)];

158: elseif endv(1)==5 159: a(n+1,1:n)=x(n)-x’; b(n)=0;

160: a(n+2,1:n)=1/2*(x(n)-x’).^2; 161: a(n+2,n+2)=x(n)-x(1);

162: else

163: error(... 164: ’Invalid value of endv in function spcof’)

165: end

166: if endv(1)==1 & n<4, c=pinv(a)*b; 167: else, c=a\b; end

168:

169: %==============================================

170:

171: function a=powermat(x,X,p) 172: %

173: % a=powermat(x,X,p)

174: % ~~~~~~~~~~~~~~~~ 175: % This function evaluates various powers of a

176: % matrix used in cubic spline interpolation. 177: %

178: % x,X - arbitrary vectors of length n and N

179: % a - an n by M matrix of elements such that 180: % a(i,j)=(x(i)>X(j))*abs(x(i)-X(j))^p

181:

© 2003 by CRC Press LLC

182: x=x(:); n=length(x); X=X(:)’; N=length(X);

183: a=x(:,ones(1,N))-X(ones(n,1),:); a=a.*(a>0);

184: switch p, case 0, a=sign(a); case 1, return; 185: case 2, a=a.*a; case 3; a=a.*a.*a;

186: case 4, a=a.*a; a=a.*a; otherwise, a=a.^p; end

4.2.3 Generalizing the Intrinsic Spline Function in MATLAB

The intrinsic MATLAB function spline employs an auxiliary function unmk to

create the piecewise polynomial definitions defining the spline. The polynomials

can be differentiated or integrated, and then functions mkpp and ppval can be used

to evaluate results. We have employed the ideas from those routines to develop func￾tions splineg and splincof extending the minimal spline capabilities of MATLAB.

The function splincof(xd,yd,endc) computes arrays b and c usable by mkpp and pp￾val. The data vector endc defines the first four types of end conditions discussed

above. The function splineg(xd,yd,x,deriv,endc,b,c) handles the same kind of data

as function spterp. Sometimes arrays b and c may have been created from a previ￾ous call to splineg or spterp. Whenever these are passed through the call list, they

are used by splineg without recomputation. Readers wanting more details on spline

concepts should consult de Boorís book [7].

Two examples illustrating spline interpolation are presented next. In the first pro￾gram called, sinetrp, a series of equally spaced points between 0 and 2π is used to

approximate y = sin(x) which satisfies

y

(x) = cos(x) , y(x) = − sin(x) ,

x

0

y(x)dx = 1 − cos(x).

The approximations for the function, derivatives, and the integral are evaluated using

splineg. Results shown in Figure 4.1 are quite satisfactory, except for points outside

the data interval [0, 2π].

© 2003 by CRC Press LLC

−0.5 0 0.5 1 1.5 2 2.5 −1.5

−1

−0.5

0

0.5

1

1.5

2

Spline Differentiation and Integration of sin(x)

x / pi

function values

y=sin(x)

data

yí(x)

yíí(x)

∫ y(x) dx

Figure 4.1: Spline Differentiation and Integration of sin(x)

© 2003 by CRC Press LLC