Attempt to code a pseudo-code. I am almost there.

I am back here hoping that someone with eagles eyes can see if I am missing a index, because I cannot find where I am missing something. I am trying to implement the following pseudo-code

Until now, my attempt to do it is:

clear all; clc;
% Parameters
n = 9; h = 1 / (n + 1);
f = @(x) 2 * (pi^2) * sin(pi * x);p = @(x) 1;q = @(x) pi^2;

% Define the function S by
S = @(x) ( (x <= -2) ) .* 0 + …
( (x > -2) & (x <= -1) ) .* (1/4 * (2 + x).^3) + …
( (x > -1) & (x <= 0) ) .* (1/4 *((2 + x).^3 – 4 * (1 + x).^3)) + …
( (x > 0) & (x <= 1) ) .* (1/4 *((2 – x).^3 – 4 * (1 – x).^3)) + …
( (x > 1) & (x <= 2) ) .* (1/4 * (2 – x).^3) + …
( (x > 2) ) .* 0;

% Define the derivative of S by
dS = @(x) ( (x <= -2) ) .* 0 + …
( (x > -2) & (x <= -1) ) .* (3/4 * (2 + x).^2) + …
( (x > -1) & (x <= 0) ) .* (3/4 *((2 + x).^2 – 4 * (1 + x).^2)) + …
( (x > 0) & (x <= 1) ) .* ((-3)/4 *((2 – x).^2 – 4 * (1 – x).^2)) + …
( (x > 1) & (x <= 2) ) .* ((-3)/4 * (2 – x).^2) + …
( (x > 2) ) .* 0;

% Difine the cubic basis splines functions by
phi = cell(n + 2, 1);
phi{1} = @(x) S( x / h ) – 4 * S( (x + h) / h );
phi{2} = @(x) S((x – h) / h) – S((x + h) / h);
for i = 3:n
phi{i} = @(x) S((x – (i-1)*h) / h);
end
phi{n + 1} = @(x) S((x – n*h) / h) – S((x – (n + 2)*h) / h);
phi{n + 2} = @(x) S((x – (n + 1)*h) / h) – 4 * S((x – (n + 2)*h) / h);

% Difine the derivatives of cubic basis splines functions by
dphi = cell(n + 2, 1);
dphi{1} = @(x) (1 / h) * dS( x / h ) – (4 / h) * dS( (x + h) / h );
dphi{2} = @(x) (1 / h) * dS((x – h) / h) – (1 / h) * dS( (x + h) / h );
for i = 3:n
dphi{i} = @(x) (1 / h) * dS((x – (i-1)*h) / h);
end
dphi{n + 1} = @(x) (1 / h) * dS((x – n*h) / h) – (1 / h) * dS((x – (n + 2)*h) / h);
dphi{n + 2} = @(x) (1 / h) * dS((x – (n + 1)*h) / h) – (4 / h) * dS((x – (n + 2)*h) / h);

% inicialize the matrix of the linear sistem and its constant vector
A = zeros(n + 2, n + 2);
b = zeros(n + 2, 1);

% Start to build A and b
for i = 1:n+2
for j= i:min(i+3,n+2)
L = max((j-3)*h, 0); % in my mind, L and U should be ajusted to this, since i cannot start from i=0
U = min((i+1)*h, 1);
integrando_aij = @(x) (p(x) * dphi{i}(x) * dphi{j}(x)) + (q(x) * phi{i}(x) * phi{j}(x));
aij = gauss_quad(integrando_aij, L, U); % the gauss_quad function is at the end of this post. It seems to work fine to me.
A(i,j) = aij;
if i~=j
A(j,i) = aij;
end
Ld = max((i-3)*h, 0); % again, the same argument as before
Ud = min((i+1)*h, 1);
integrando_di = @(x) f(x) * phi{i}(x);
di = gauss_quad(integrando_di, Ld, Ud);
b(i) = di;
end
end

c = A b;

disp(‘Coefficients:’);
for i = 0:n+1
fprintf(‘c%d = %.6fn’, i, c(i+1));
end

%%%%%%%%%% PLOTTING GRAPHS %%%%%%%%%%%%%

% Here I just want to look to the shape of the functions that I have
% defined, just to make sure I am using rights tools

% Number of pts to plot
num_points = 1000;
x_values = linspace(0, 1, num_points);

figure (1);
hold on;

% Plot each function phi in [0, 1]
for i = 1:n+2
y_values = arrayfun(phi{i}, x_values);
plot(x_values, y_values, ‘DisplayName’, sprintf(‘\phi_%d(x)’, i-1));
end
xlabel(‘x’);
ylabel(‘phi_i(x)’);
title(‘Funções base phi_i(x) entre [0, 1]’);
legend show;
hold off;

figure (2)
% Plot an specific function, here I chose phi{11}
y_values_phi = arrayfun(phi{11}, x_values);
plot(x_values, y_values_phi);
xlabel(‘x’);
ylabel(‘phi(x)’);
title(‘Função base phi(x) entre [0, 1]’);
legend show;

figure (3)
% Plot S and its derivative
x_values_s = linspace(-3, 3, num_points);
y_values_s = arrayfun(S, x_values_s);
plot(x_values_s, y_values_s);
xlabel(‘x’);
ylabel(‘y’);
hold on
y_values_ds = arrayfun(dS, x_values_s);
plot(x_values_s, y_values_ds);
legend(‘S(x)’, ‘dS(x)’)
hold off

In the book that I am following the basis functions shoul look like

The last section of my code, suggest that I am using them right. Finally, here is the gauss_quad function that I call to solve integrals
function result = gauss_quad(f, a, b)
% 3-point Gauss quadrature
x = [-sqrt(3/5), 0, sqrt(3/5)];
w = [5/9, 8/9, 5/9];

% Transform the integration limits
t = 0.5 * ((b – a) * x + (b + a));
ft = arrayfun(f, t);

result = 0.5 * (b – a) * sum(w .* ft);
end
The coefficients should be
.
I think I have done a good work so far, and I need someone to give me a direction.I am back here hoping that someone with eagles eyes can see if I am missing a index, because I cannot find where I am missing something. I am trying to implement the following pseudo-code

Until now, my attempt to do it is:

clear all; clc;
% Parameters
n = 9; h = 1 / (n + 1);
f = @(x) 2 * (pi^2) * sin(pi * x);p = @(x) 1;q = @(x) pi^2;

% inicialize the matrix of the linear sistem and its constant vector
A = zeros(n + 2, n + 2);
b = zeros(n + 2, 1);

c = A b;

disp(‘Coefficients:’);
for i = 0:n+1
fprintf(‘c%d = %.6fn’, i, c(i+1));
end

%%%%%%%%%% PLOTTING GRAPHS %%%%%%%%%%%%%

% Here I just want to look to the shape of the functions that I have
% defined, just to make sure I am using rights tools

% Number of pts to plot
num_points = 1000;
x_values = linspace(0, 1, num_points);

figure (1);
hold on;

In the book that I am following the basis functions shoul look like

% Transform the integration limits
t = 0.5 * ((b – a) * x + (b + a));
ft = arrayfun(f, t);

result = 0.5 * (b – a) * sum(w .* ft);
end
The coefficients should be
.
I think I have done a good work so far, and I need someone to give me a direction. I am back here hoping that someone with eagles eyes can see if I am missing a index, because I cannot find where I am missing something. I am trying to implement the following pseudo-code

Until now, my attempt to do it is:

clear all; clc;
% Parameters
n = 9; h = 1 / (n + 1);
f = @(x) 2 * (pi^2) * sin(pi * x);p = @(x) 1;q = @(x) pi^2;

% inicialize the matrix of the linear sistem and its constant vector
A = zeros(n + 2, n + 2);
b = zeros(n + 2, 1);

c = A b;

disp(‘Coefficients:’);
for i = 0:n+1
fprintf(‘c%d = %.6fn’, i, c(i+1));
end

%%%%%%%%%% PLOTTING GRAPHS %%%%%%%%%%%%%

% Here I just want to look to the shape of the functions that I have
% defined, just to make sure I am using rights tools

% Number of pts to plot
num_points = 1000;
x_values = linspace(0, 1, num_points);

figure (1);
hold on;

In the book that I am following the basis functions shoul look like

% Transform the integration limits
t = 0.5 * ((b – a) * x + (b + a));
ft = arrayfun(f, t);

result = 0.5 * (b – a) * sum(w .* ft);
end
The coefficients should be
.
I think I have done a good work so far, and I need someone to give me a direction. code, homework, algorithm, integration, index MATLAB Answers — New Questions

Cart

Cart