% CUBIC SPLINE RAYLEIGH-RITZ ALGORITHM 11.6

%

% 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)=0

%

% 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))**3

% 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+1

% since PHI(I) is nonzero only between X(I-2) and X(I+2)

% unless I = 1, 2, N+1 or N+2

% (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])^3

% 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');

docsity.com

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,' pi*pi \n');

fprintf(1,' 2*pi*pi*sin(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;

docsity.com

N2 = N+2;

N3 = N+3;

% Initialize matrix A at zero, note that A(I, N+3) =B(I)

A = zeros(N2,N3);

C = zeros(1,N+2);

X = zeros(1,N+2);

CO = zeros(N+2,4,4);

DCO = zeros(N+2,4,3);

AF = zeros(1,N+1);

BF = zeros(1,N+1);

CF = zeros(1,N+1);

DF = zeros(1,N+1);

AP = zeros(1,N+1);

BP = zeros(1,N+1);

CP = zeros(1,N+1);

DP = zeros(1,N+1);

AQ = zeros(1,N+1);

BQ = zeros(1,N+1);

CQ = zeros(1,N+1);

DQ = zeros(1,N+1);

AA1 = zeros(1,N+2);

BB1 = zeros(1,N+2);

CC1 = zeros(1,N+2);

DD1 = zeros(1,N+2);

XA = zeros(1,N+2);

XL1 = zeros(1,N+2);

XU1 = zeros(1,N+2);

XZ = zeros(1,N+2);

% STEP 2

% X(1)=0,...,X(I) = (I-1)*H,...,X(N+1) = 1-H, X(N+2) = 1

for I = 1 : N2

X(I) = (I-1)*H;

end;

% STEPS 3 and 4 are implemented in what follows. Initialize coefficients

% CO(I,J,K), DCO(I,J,K) */

for I = 1 : N2

for J = 1 : 4

% JJ corresponds the coefficients of phi and phi' to the proper interval

% involving J */

JJ = I+J-3;

CO(I,J,1) = 0;

CO(I,J,2) = 0;

CO(I,J,3) = 0;

CO(I,J,4) = 0;

E = I-1;

OK = TRUE;

if JJ < I-2 | JJ >= I+2

OK = FALSE;

end;

if I == 1 & JJ < I

OK = FALSE;

end;

if I == 2 & JJ < I-1

OK = FALSE;

docsity.com

end;

if I == N+1 & JJ > N+1

OK = FALSE;

end;

if I == N+2 & JJ >= N+2

OK = FALSE;

end;

if OK == TRUE

if JJ <= I-2

CO(I,J,1) = (((-E+6)*E-12)*E+8)/24;

CO(I,J,2) = ((E-4)*E+4)/(8*H);

CO(I,J,3) = (-E+2)/(8*H^2);

CO(I,J,4) = 1/(24*H^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)/(8*H);

CO(I,J,3) = (E+2)/(8*H^2);

CO(I,J,4) = -1/(24*H^3);

OK = FALSE;

else

if JJ > I-1

CO(I,J,1) = ((-3*E-6)*E*E+4)/24;

CO(I,J,2) = (3*E+4)*E/(8*H);

CO(I,J,3) = (-3*E-2)/(8*H^2);

CO(I,J,4) = 1/(8*H^3);

if I ~= 1 & I ~= N+1

OK = FALSE;

end;

else

CO(I,J,1) = ((3*E-6)*E*E+4)/24;

CO(I,J,2) = (-3*E+4)*E/(8*H);

CO(I,J,3) = (3*E-2)/(8*H^2);

CO(I,J,4) = -1/(8*H^3);

if I ~= 2 & I ~= N+2

OK = FALSE;

end;

end;

end;

end;

end;

if OK == TRUE

if I <= 2

AA = 1/24;

BB = -1/(8*H);

CC = 1/(8*H^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

docsity.com

##### Document information

Uploaded by:
saripella

Views: 2687

Downloads :
11

University:
Biyani Girls College

Subject:
Numerical Analysis

Upload date:
31/07/2012

very goooood