# 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  -  4  of  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');
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