Download Boundary Value Problem-Numerical Analysis-MATLAB Code and more Exercises Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
% CUBIC SPLINE RAYLEIGH-RITZ ALGORITHM 11.
% To approximate the solution to the boundary-value problem % % -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1, Y(0)=Y(1)= % % With a sum of cubic splines: % % INPUT: Integer n % % OUTPUT: Coefficients C(0),...,C(n+1) of the basis functions % % GENERAL OUTLINE % % 1. Nodes labelled X(I)=(I-1)H, 1 <= I <= N+2, where % H=1/(N+1) so that zero subscripts are avoided % 2. The functions PHI(I) and PHI'(I) are shifted so that % PHI(1) and PHI'(1) are centered at X(1), PHI(2) and PHI'(2) % are centered at X(2),... , PHI(N+2) and % PHI'(N+2) are centered at (X(N+2)---for example, % PHI(3) = S((X-X(3))/H) % = S(X/H + 2) % 3. The functions PHI(I) are represented in terms of their % coefficients in the following way: % (PHI(I))(X) = CO(I,K,1) + CO(I,K,2)(X-X(J)) + % CO(I,K,3)(X-X(J))2 + CO(I,K,4)(X-X(J))** % for X(J) <= X <= X(J+1) where % K=1 IF J=I-2, K=2 IF J=I-1, K=3 IF J=I, K=4 IF J=I+ % since PHI(I) is nonzero only between X(I-2) and X(I+2) % unless I = 1, 2, N+1 or N+ % (see subroutine PHICO) % 4. The derivative of PHI(I) denoted PHI'(I) is represented % as in 3. By its coefficients DCO(I,K,L), L = 1, 2, 3 % (See subroutine DPHICO). % 5. The functions P,Q and F are represented by their cubic % spline interpolants using clamped boundary conditions % (see Algorithm 3.5). Thus, for X(I) <= X <= X(I+1) we % use AF(I)+BF(I)(X-X[I])+CF(I)(X-X[I])^2+DF(I)*(X-X[I])^ % to represent F(X). Similarly, AP,BP,CP,DP are used for P % and AQ,BQ,CQ,DQ are used for Q. (see subroutine COEF). % 6. The integrands in STEPS 6 and 9 are replaced by products % of cubic polynomial approximations on each subinterval of % length H and the integrals of the resulting polynomials % are computed exactly. (see subroutine XINT). % % syms('s','S','SS','FPL','FPR','PPL','PPR','QPL','QPR','OK'); syms('N','FLAG','NAME','OUP','H','N1','N2','N3','A','C','X'); syms('CO','DCO','AF','BF','CF','DF','AP','BP','CP','DP'); syms('AQ','BQ','CQ','DQ','AA','BB','CC','DD','AA1','BB1'); syms('CC1','DD1','XU1','XL1','XA','XL','XZ','XU','I','J','K'); syms('JJ','J0','J1','J2','JJ1','JJ2','KK','E','A1','B1','C1','D1'); syms('A2','B2','C2','D2','A3','B3','C3','D3','A4','B4','C4','D4'); syms('ZZ1','ZZ2','K2','K3');
TRUE = 1;
FALSE = 0;
fprintf(1,'This is the Cubic Spline Rayleigh-Ritz Method.\n'); OK = FALSE; fprintf(1,'Input functions P(X), Q(X), and F(X) on separate lines.\n'); fprintf(1,'For example: 1 \n'); fprintf(1,' pipi \n'); fprintf(1,' 2pipisin(pi*x) \n'); s = input(' ','s'); P = inline(s,'x'); s = input(' ','s'); Q = inline(s,'x'); s = input(' ','s'); F = inline(s,'x'); fprintf(1,'Input derivative of F evaluated at 0 \n'); FPL = input(' '); fprintf(1,'Input derivative of F evaluated at 1 \n'); FPR = input(' '); fprintf(1,'Input derivative of Q evaluated at 0 \n'); QPL = input(' '); fprintf(1,'Input derivative of Q evaluated at 1 \n'); QPR = input(' '); fprintf(1,'Input derivative of P evaluated at 0 \n'); PPL = input(' '); fprintf(1,'Input derivative of P evaluated at 1 \n'); PPR = input(' '); while OK == FALSE fprintf(1,'Input positive integer n, where x(0) = 0, '); fprintf(1,'..., x(n+1) = 1.\n'); N = input(' '); if N <= 0 fprintf(1,'Number must be a positive integer.\n'); else OK = TRUE; end; end; if OK == TRUE fprintf(1,'Choice of output method:\n'); fprintf(1,'1. Output to screen\n'); fprintf(1,'2. Output to text File\n'); fprintf(1,'Please enter 1 or 2.\n'); FLAG = input(' '); if FLAG == 2 fprintf(1,'Input the file name in the form - drive:\name.ext\n'); fprintf(1,'for example A:\OUTPUT.DTA\n'); NAME = input(' ','s'); OUP = fopen(NAME,'wt'); else OUP = 1; end; fprintf(OUP, 'CUBIC SPLINE RAYLEIGH-RITZ METHOD\n\n'); % STEP 1 H = 1/(N+1); N1 = N+1;
end; if I == N+1 & JJ > N+ OK = FALSE; end; if I == N+2 & JJ >= N+ OK = FALSE; end; if OK == TRUE if JJ <= I- CO(I,J,1) = (((-E+6)E-12)E+8)/24; CO(I,J,2) = ((E-4)E+4)/(8H); CO(I,J,3) = (-E+2)/(8H^2); CO(I,J,4) = 1/(24H^3); OK = FALSE; else if JJ > I CO(I,J,1) = (((E+6)E+12)E+8)/24; CO(I,J,2) = ((-E-4)E-4)/(8H); CO(I,J,3) = (E+2)/(8H^2); CO(I,J,4) = -1/(24H^3); OK = FALSE; else if JJ > I- CO(I,J,1) = ((-3E-6)EE+4)/24; CO(I,J,2) = (3E+4)E/(8H); CO(I,J,3) = (-3E-2)/(8H^2); CO(I,J,4) = 1/(8H^3); if I ~= 1 & I ~= N+ OK = FALSE; end; else CO(I,J,1) = ((3E-6)EE+4)/24; CO(I,J,2) = (-3E+4)E/(8H); CO(I,J,3) = (3E-2)/(8H^2); CO(I,J,4) = -1/(8H^3); if I ~= 2 & I ~= N+ OK = FALSE; end; end; end; end; end; if OK == TRUE if I <= 2 AA = 1/24; BB = -1/(8H); CC = 1/(8H^2); DD = -1/(24*H^3); if I == 2 CO(I,J,1) = CO(I,J,1)-AA; CO(I,J,2) = CO(I,J,2)-BB; CO(I,J,3) = CO(I,J,3)-CC; CO(I,J,4) = CO(I,J,4)-DD; else
CO(I,J,1) = CO(I,J,1)-4*AA;
CO(I,J,2) = CO(I,J,2)-4*BB;
CO(I,J,3) = CO(I,J,3)-4*CC;
CO(I,J,4) = CO(I,J,4)-4*DD;
end; else EE = N+2; AA = (((-EE+6)EE-12)EE+8)/24; BB = ((EE-4)EE+4)/(8H); CC = (-EE+2)/(8H^2); DD = 1/(24H^3); if I == N+ CO(I,J,1) = CO(I,J,1)-AA; CO(I,J,2) = CO(I,J,2)-BB; CO(I,J,3) = CO(I,J,3)-CC; CO(I,J,4) = CO(I,J,4)-DD; else CO(I,J,1) = CO(I,J,1)-4AA; CO(I,J,2) = CO(I,J,2)-4BB; CO(I,J,3) = CO(I,J,3)-4CC; CO(I,J,4) = CO(I,J,4)-4DD; end; end; end; DCO(I,J,1) = 0; DCO(I,J,2) = 0; DCO(I,J,3) = 0; E = I-1; OK = TRUE; if JJ < I-2 | JJ >= I+ OK = FALSE; end; if I == 1 & JJ < I OK = FALSE; end; if I == 2 & JJ < I- OK = FALSE; end; if I == N+1 & JJ > N+ OK = FALSE; end; if I == N+2 & JJ >= N+ OK = FALSE; end; if OK == TRUE if JJ <= I- DCO(I,J,1) = ((E-4)E+4)/(8H); DCO(I,J,2) = (-E+2)/(4H^2); DCO(I,J,3) = 1/(8H^3); OK = FALSE; else if JJ > I DCO(I,J,1) = ((-E-4)E-4)/(8H);
% Output the basis functions. / fprintf(OUP, 'Basis Function: A + BX + CX2 + DX3\n\n'); fprintf(OUP, ' A B C D\n\n'); for I = 1 : N fprintf(OUP, 'phi( %d )\n\n', I); for J = 1 : 4 if I ~= 1 | (J ~= 1 & J ~= 2) if I ~= 2 | J ~= 2 if I ~= N1 | J ~= 4 if I ~= N2 | (J ~= 3 & J ~= 4) JJ1 = I+J-3; JJ2 = I+J-2; fprintf(OUP, 'On (X( %d ), X( %d )) ', JJ1, JJ2); for K = 1 : 4 fprintf(OUP, ' %12.8f ', CO(I,J,K)); end; fprintf(OUP, '\n'); end; end; end; end; end; end; % Obtain coefficients for F, P, Q for I = 1 : N AA1(I) = F(X(I)); end; XA(1) = 3.0(AA1(2)-AA1(1))/H-3.0FPL; XA(N2) = 3.0FPR-3.0(AA1(N2)-AA1(N2-1))/H; XL1(1) = 2.0H; XU1(1) = 0.5; XZ(1) = XA(1)/XL1(1); for I = 2 : N XA(I) = 3.0(AA1(I+1)-2.0AA1(I)+AA1(I-1))/H; XL1(I) = H(4.0-XU1(I-1)); XU1(I) = H/XL1(I); XZ(I) = (XA(I)-HXZ(I-1))/XL1(I); end; XL1(N2) = H(2.0-XU1(N2-1)); XZ(N2) = (XA(N2)-HXZ(N2-1))/XL1(N2); CC1(N2) = XZ(N2); for I = 1 : N J = N2-I; CC1(J) = XZ(J)-XU1(J)CC1(J+1); BB1(J) = (AA1(J+1)-AA1(J))/H-H(CC1(J+1)+2.0CC1(J))/3.0; DD1(J) = (CC1(J+1)-CC1(J))/(3.0H); end; for I = 1 : N AF(I) = ((-DD1(I)X(I)+CC1(I))X(I)-BB1(I))X(I)+AA1(I); BF(I) = (3.0DD1(I)X(I)-2.0CC1(I))X(I)+BB1(I); CF(I) = CC1(I)-3.0DD1(I)X(I); DF(I) = DD1(I); end; for I = 1 : N
AA1(I) = P(X(I));
end; XA(1) = 3.0(AA1(2)-AA1(1))/H-3.0PPL; XA(N2) = 3.0PPR-3.0(AA1(N2)-AA1(N2-1))/H; XL1(1) = 2.0H; XU1(1) = 0.5; XZ(1) = XA(1)/XL1(1); for I = 2 : N XA(I) = 3.0(AA1(I+1)-2.0AA1(I)+AA1(I-1))/H; XL1(I) = H(4.0-XU1(I-1)); XU1(I) = H/XL1(I); XZ(I) = (XA(I)-HXZ(I-1))/XL1(I); end; XL1(N2) = H(2.0-XU1(N2-1)); XZ(N2) = (XA(N2)-HXZ(N2-1))/XL1(N2); CC1(N2) = XZ(N2); for I = 1 : N J = N2-I; CC1(J) = XZ(J)-XU1(J)CC1(J+1); BB1(J) = (AA1(J+1)-AA1(J))/H -H(CC1(J+1)+2.0CC1(J))/3.0; DD1(J) = (CC1(J+1)-CC1(J))/(3.0H); end; for I = 1 : N AP(I) = ((-DD1(I)X(I)+CC1(I))X(I)-BB1(I))X(I)+AA1(I); BP(I) = (3.0DD1(I)X(I)-2.0CC1(I))X(I)+BB1(I); CP(I) = CC1(I)-3.0DD1(I)X(I); DP(I) = DD1(I); end; for I = 1 : N AA1(I) = Q(X(I)); end; XA(1) = 3.0(AA1(2)-AA1(1))/H-3.0QPL; XA(N2) = 3.0QPR-3.0(AA1(N2)-AA1(N2-1))/H; XL1(1) = 2.0H; XU1(1) = 0.5; XZ(1) = XA(1)/XL1(1); for I = 2 : N XA(I) = 3.0(AA1(I+1)-2.0AA1(I)+AA1(I-1))/H; XL1(I) = H(4.0-XU1(I-1)); XU1(I) = H/XL1(I); XZ(I) = (XA(I)-HXZ(I-1))/XL1(I); end; XL1(N2) = H(2.0-XU1(N2-1)); XZ(N2) = (XA(N2)-HXZ(N2-1))/XL1(N2); CC1(N2) = XZ(N2); for I = 1 : N J = N2-I; CC1(J) = XZ(J)-XU1(J)CC1(J+1); BB1(J) = (AA1(J+1)-AA1(J))/H -H(CC1(J+1)+2.0CC1(J))/3.0; DD1(J) = (CC1(J+1)-CC1(J))/(3.0H); end; for I = 1 : N AQ(I) = ((-DD1(I)X(I)+CC1(I))X(I)-BB1(I))X(I)+AA1(I); BQ(I) = (3.0DD1(I)X(I)-2.0CC1(I))X(I)+BB1(I);
A4 = CO(J,K,1);
B4 = CO(J,K,2);
C4 = CO(J,K,3);
D4 = CO(J,K,4);
ZZ1 = XINT(XU,XL,AP(JJ),BP(JJ),CP(JJ),DP(JJ),A1,B1,C1,D1,A3,B3,C3,D3);
ZZ2 = XINT(XU,XL,AQ(JJ),BQ(JJ),CQ(JJ),DQ(JJ),A2,B2,C2,D2,A4,B4,C4,D4);
A(I,J) = A(I,J) + ZZ1 + ZZ2;
end; A(J,I) = A(I,J); end; end; end; % STEP 10 for I = 1 : N for J = I+1 : N CC = A(J,I)/A(I,I); for K = I+1 : N A(J,K) = A(J,K)-CCA(I,K); end; A(J,I) = 0; end; end; C(N2) = A(N2,N3)/A(N2,N2); for I = 1 : N J = N1-I+1; C(J) = A(J,N3); for KK = J+1 : N C(J) = C(J)-A(J,KK)C(KK); end; C(J) = C(J)/A(J,J); end; % STEP 11 % output coefficients fprintf(OUP, '\nCoefficients: c(1), c(2), ... , c(n+1)\n\n'); for I = 1 : N fprintf(OUP, ' %12.6e \n', C(I)); end; fprintf(OUP, '\n'); % compute and output value of the approximation at the nodes fprintf(OUP, 'The approximation evaluated at the nodes:\n\n'); fprintf(OUP, ' Node Value\n\n'); for I = 1 : N S = 0; for J = 1 : N J0 = max(J-2,1); J1 = min(J+2,N+2); SS = 0; if I < J0 | I >= J S = S + C(J)SS; else K = INTE(J,I); SS = ((CO(J,K,4)X(I)+CO(J,K,3))X(I)+CO(J,K,2))X(I)+CO(J,K,1); S = S + C(J)*SS; end;
end; fprintf(OUP, '%12.8f %12.8f\n', X(I), S); end; end; if OUP ~= 1 fclose(OUP); fprintf(1,'Output file %s created successfully \n',NAME); end;
