Boundary Value Problem-Numerical Analysis-MATLAB Code, Exercises for Numerical Analysis. Biyani Girls College

Numerical Analysis

Description: This is solution to one of problems in Numerical Analysis. This is matlab code. Its helpful to students of Computer Science, Electrical and Mechanical Engineering. This code also help to understand algorithm and logic behind the problem. This code includes: Boundary, Value, Problem, Cubic, Spline, Rayleigh, Ritz, Algorithm, Coefficients, Basis, Function, Derivative, Interpolants
Showing pages  1  -  2  of  11
The preview of this document ends here! Please or to read the full document or to download it.
Document information
Embed this document:
aaxcxaa - Amity University

very goooood

20/02/13 02:11

% 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

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)/(8*H);

CC = (-EE+2)/(8*H^2);

DD = 1/(24*H^3);

if I == N+1

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;

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+2

OK = FALSE;

end;

if I == 1 & JJ < I

OK = FALSE;

end;

if I == 2 & JJ < I-1

OK = FALSE;

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

DCO(I,J,1) = ((E-4)*E+4)/(8*H);

DCO(I,J,2) = (-E+2)/(4*H^2);

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

OK = FALSE;

else

if JJ > I

DCO(I,J,1) = ((-E-4)*E-4)/(8*H);

docsity.com

DCO(I,J,2) = (E+2)/(4*H^2);

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

OK = FALSE;

else

if JJ > I-1

DCO(I,J,1) = (3*E+4)*E/(8*H);

DCO(I,J,2) = (-3.0*E-2.0)/(4.0*H^2);

DCO(I,J,3) = 3/(8*H^3);

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

OK = FALSE;

end;

else

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

DCO(I,J,2) = (3*E-2)/(4*H^2);

DCO(I,J,3) = -3/(8*H^3);

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

OK = FALSE;

end;

end;

end;

end;

end;

if OK == TRUE

if I <= 2

AA = -1/(8*H);

BB = 1/(4*H^2);

CC = -1/(8*H^3);

if I == 2

DCO(I,J,1) = DCO(I,J,1)-AA;

DCO(I,J,2) = DCO(I,J,2)-BB;

DCO(I,J,3) = DCO(I,J,3)-CC;

else

DCO(I,J,1) = DCO(I,J,1)-4*AA;

DCO(I,J,2) = DCO(I,J,2)-4*BB;

DCO(I,J,3) = DCO(I,J,3)-4*CC;

end;

else

EE = N+2;

AA = ((EE-4)*EE+4)/(8*H);

BB = (-EE+2)/(4*H^2);

CC = 1/(8*H^3);

if I == N+1

DCO(I,J,1) = DCO(I,J,1)-AA;

DCO(I,J,2) = DCO(I,J,2)-BB;

DCO(I,J,3) = DCO(I,J,3)-CC;

else

DCO(I,J,1) = DCO(I,J,1)-4*AA;

DCO(I,J,2) = DCO(I,J,2)-4*BB;

DCO(I,J,3) = DCO(I,J,3)-4*CC;

end;

end;

end;

end;

end;

docsity.com

% Output the basis functions. */

fprintf(OUP, 'Basis Function: A + B*X + C*X**2 + D*X**3\n\n');

fprintf(OUP, ' A B C D\n\n');

for I = 1 : N2

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 : N2

AA1(I) = F(X(I));

end;

XA(1) = 3.0*(AA1(2)-AA1(1))/H-3.0*FPL;

XA(N2) = 3.0*FPR-3.0*(AA1(N2)-AA1(N2-1))/H;

XL1(1) = 2.0*H;

XU1(1) = 0.5;

XZ(1) = XA(1)/XL1(1);

for I = 2 : N1

XA(I) = 3.0*(AA1(I+1)-2.0*AA1(I)+AA1(I-1))/H;

XL1(I) = H*(4.0-XU1(I-1));

XU1(I) = H/XL1(I);

XZ(I) = (XA(I)-H*XZ(I-1))/XL1(I);

end;

XL1(N2) = H*(2.0-XU1(N2-1));

XZ(N2) = (XA(N2)-H*XZ(N2-1))/XL1(N2);

CC1(N2) = XZ(N2);

for I = 1 : N1

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.0*CC1(J))/3.0;

DD1(J) = (CC1(J+1)-CC1(J))/(3.0*H);

end;

for I = 1 : N1

AF(I) = ((-DD1(I)*X(I)+CC1(I))*X(I)-BB1(I))*X(I)+AA1(I);

BF(I) = (3.0*DD1(I)*X(I)-2.0*CC1(I))*X(I)+BB1(I);

CF(I) = CC1(I)-3.0*DD1(I)*X(I);

DF(I) = DD1(I);

end;

for I = 1 : N2

docsity.com

AA1(I) = P(X(I));

end;

XA(1) = 3.0*(AA1(2)-AA1(1))/H-3.0*PPL;

XA(N2) = 3.0*PPR-3.0*(AA1(N2)-AA1(N2-1))/H;

XL1(1) = 2.0*H;

XU1(1) = 0.5;

XZ(1) = XA(1)/XL1(1);

for I = 2 : N1

XA(I) = 3.0*(AA1(I+1)-2.0*AA1(I)+AA1(I-1))/H;

XL1(I) = H*(4.0-XU1(I-1));

XU1(I) = H/XL1(I);

XZ(I) = (XA(I)-H*XZ(I-1))/XL1(I);

end;

XL1(N2) = H*(2.0-XU1(N2-1));

XZ(N2) = (XA(N2)-H*XZ(N2-1))/XL1(N2);

CC1(N2) = XZ(N2);

for I = 1 : N1

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.0*CC1(J))/3.0;

DD1(J) = (CC1(J+1)-CC1(J))/(3.0*H);

end;

for I = 1 : N1

AP(I) = ((-DD1(I)*X(I)+CC1(I))*X(I)-BB1(I))*X(I)+AA1(I);

BP(I) = (3.0*DD1(I)*X(I)-2.0*CC1(I))*X(I)+BB1(I);

CP(I) = CC1(I)-3.0*DD1(I)*X(I);

DP(I) = DD1(I);

end;

for I = 1 : N2

AA1(I) = Q(X(I));

end;

XA(1) = 3.0*(AA1(2)-AA1(1))/H-3.0*QPL;

XA(N2) = 3.0*QPR-3.0*(AA1(N2)-AA1(N2-1))/H;

XL1(1) = 2.0*H;

XU1(1) = 0.5;

XZ(1) = XA(1)/XL1(1);

for I = 2 : N1

XA(I) = 3.0*(AA1(I+1)-2.0*AA1(I)+AA1(I-1))/H;

XL1(I) = H*(4.0-XU1(I-1));

XU1(I) = H/XL1(I);

XZ(I) = (XA(I)-H*XZ(I-1))/XL1(I);

end;

XL1(N2) = H*(2.0-XU1(N2-1));

XZ(N2) = (XA(N2)-H*XZ(N2-1))/XL1(N2);

CC1(N2) = XZ(N2);

for I = 1 : N1

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.0*CC1(J))/3.0;

DD1(J) = (CC1(J+1)-CC1(J))/(3.0*H);

end;

for I = 1 : N1

AQ(I) = ((-DD1(I)*X(I)+CC1(I))*X(I)-BB1(I))*X(I)+AA1(I);

BQ(I) = (3.0*DD1(I)*X(I)-2.0*CC1(I))*X(I)+BB1(I);

docsity.com

CQ(I) = CC1(I)-3.0*DD1(I)*X(I);

DQ(I) = DD1(I);

end;

% STEPS 5-9 are implemented in what follows

for I = 1 : N2

% indices for limits of integration for A(I,J) and B(I)

J1 = min(I+2,N+2);

J0 = max(I-2,1);

J2 = J1-1;

% integrate over each subinterval where phi(I) is nonzero

for JJ = J0 : J2

% limits of integration for each call

XU = X(JJ+1);

XL = X(JJ);

% coefficients of bases

K = INTE(I,JJ);

A1 = DCO(I,K,1);

B1 = DCO(I,K,2);

C1 = DCO(I,K,3);

D1 = 0;

A2 = CO(I,K,1);

B2 = CO(I,K,2);

C2 = CO(I,K,3);

D2 = CO(I,K,4);

% call subprogram for integrations

ZZ1 = XINT(XU,XL,AP(JJ),BP(JJ),CP(JJ),DP(JJ),A1,B1,C1,D1,A1,B1,C1,D1);

ZZ2 = XINT(XU,XL,AQ(JJ),BQ(JJ),CQ(JJ),DQ(JJ),A2,B2,C2,D2,A2,B2,C2,D2);

A(I,I) = A(I,I) + ZZ1 + ZZ2;

ZZ1 = XINT(XU,XL,AF(JJ),BF(JJ),CF(JJ),DF(JJ),A2,B2,C2,D2,1,0,0,0);

A(I,N+3) = A(I,N+3) + ZZ1;

end;

% compute A(I,J) for J = I+1,...,Min(I+3,N+2)

K3 = I+1;

if K3 <= N2

K2 = min(I+3,N+2);

for J = K3 : K2

J0 = max(J-2,1);

for JJ = J0 : J2

XU = X(JJ+1);

XL = X(JJ);

K = INTE(I,JJ);

A1 = DCO(I,K,1);

B1 = DCO(I,K,2);

C1 = DCO(I,K,3);

D1 = 0;

A2 = CO(I,K,1);

B2 = CO(I,K,2);

C2 = CO(I,K,3);

D2 = CO(I,K,4);

K = INTE(J,JJ);

A3 = DCO(J,K,1);

B3 = DCO(J,K,2);

C3 = DCO(J,K,3);

D3 = 0;

docsity.com

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 : N1

for J = I+1 : N2

CC = A(J,I)/A(I,I);

for K = I+1 : N3

A(J,K) = A(J,K)-CC*A(I,K);

end;

A(J,I) = 0;

end;

end;

C(N2) = A(N2,N3)/A(N2,N2);

for I = 1 : N1

J = N1-I+1;

C(J) = A(J,N3);

for KK = J+1 : N2

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 : N1

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 : N2

S = 0;

for J = 1 : N2

J0 = max(J-2,1);

J1 = min(J+2,N+2);

SS = 0;

if I < J0 | I >= J1

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;

docsity.com

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;



docsity.com

Docsity is not optimized for the browser you're using. In order to have a better experience please switch to Google Chrome, Firefox, Internet Explorer 9+ or Safari! Download Google Chrome