Compact approximation stencils based on integrated flat radial basis functions.Compact local integrated RBF stencils.

1 view (last 30 days)
clearvars
clc
f=@(x)-exp(-5*x)*(9975*sin(100*x)+1000*cos(100*x));
% parameters
nx=91;
Lx=1;
dx=Lx/(nx-1);
h=dx;
x=linspace(0,Lx,nx);
%Boundary conditions
u(1,:)=0;
u(nx,:)=0;
% Multiquadric RBF
beta=20;
a=beta*h;
% Define G(x)
G=zeros(nx,nx+2);
for j=1:nx
for i=1:nx
G(i,j)=sqrt((x(i)-x(j))^2+a^2);
end
end
% Extract the first and last rows of the matrix G
G1=G(1,:);
G2=G(end,:);
G_bar=[G1;G2];
G_a=G(2:end-1,:);
%H
H=zeros(nx,nx);
for i=1:nx
for j=1:nx
H(i,j)=(((x(i)-x(j))*sqrt((x(i)-x(j)).^2+a.^2))/2)+((a.^2)/2)*...
log((x(i)-x(j))+sqrt((x(i)-x(j)).^2+a.^2));
end
end
H_a=[H,ones(nx,1),zeros(nx,1)];
% Plot H
% surf(x, x, H); % 3D plot of H
% xlabel('x');
% ylabel('x');
% zlabel('H');
% title('Plot of H');
% colorbar;
% hold on
%H_bar
H1=zeros(nx,nx);
for i=1:nx
for j=1:nx
H1(i,j)=((sqrt((x(i)-x(j)).^2+a.^2))/6)+((a.^2)/2)*(x(i)-x(j))*....
log((x(i)-x(j))+sqrt((x(i)-x(j)).^2+a.^2))-((a.^2)/2)*...
sqrt((x(i)-x(j)).^2+a.^2);
end
end
% 3D plot of H1
% surf(x, x, H1);
% xlabel('x');
% ylabel('x');
% zlabel('H1');
% title('Plot of H1');
% colorbar;
% Add the new column x and 1
H_bar=[H1,x',ones(size(H,1),1)];
H_inv=pinv(H_bar);
% Conversion matrix
C=[H_bar;G_bar];
C_inv=pinv(C);
%f
f_1=f(0);
f_nx=f(end);
D=[u;f_1;f_nx];
E=C_inv*D;
f_interior=G_a*E;
% Populate the tri-diagonal structure of A
% Construct the system matrix A
A=zeros(nx-2,nx-2);
for i=2:nx-1
A(i-1,i-1)=f_interior(i-1);
if i>2
A(i-1,i-1)=f_interior(i-1);
end
if i<nx-1
A(i-1,i-1)=f_interior(i);
end
end
% Construct the vector b
b=zeros(nx-2,1);
for i=1:nx-1
b(i)=-exp(-5*x(i))*(9975*sin(100*x(i))+1000*cos(100*x(i)));
end
%LU decompose
N=length(A);
L=zeros(N,N);
U=zeros(N,N);
for i=1:N
L(i,i)=1;
end
U(1,:)=A(1,:);
L(:,1)=A(:,1)/U(1,1);
for i=2:N
for j=i:N
U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
end
for k=i+1:N
L(k,i)=(A(k,i)-L(k,1:i-1)*U(1:i-1,i))/U(i,i);
end
end
%Solve Ly = B for y using forward substitution
y=zeros(N,1);
y(1)=b(1)/L(1,1);
for k=2:N
y(k)=(b(k)-L(k,1:k-1)*y(1:k-1))/L(k,k);
end
% Solve Ux = y for x using backward substitution
z=zeros(N,1);
z(N)=y(N)/U(N,N);
for k=N-1:-1:1
z(k)=(y(k)-U(k,k+1:N)*z(k+1:N))/U(k,k);
end
% Plot numerical solution
% Adjust the length of x and z for plotting
x_plot = x(2:end-1); % Adjust x to match the length of z
z_plot = z; % Use z as is
% Plot numerical solution
figure;
plot(x_plot, z_plot, 'r-', 'LineWidth', 1.5);
hold on;
% Define the exact solution function
uexact=@(x)sin(100*x).*exp(-5*x);
% Evaluate the exact solution
u_exact_values=uexact(x);
% Plot the exact solution
figure;
plot(x,u_exact_values,'b','LineWidth',1.5);
xlabel('x');
ylabel('u(x)');
title('Exact Solution');
% Compute the RMS error
u_interior_exact = uexact(x(2:end-1));
RMS_error = sqrt(sum((z - u_interior_exact).^2)/nx);
disp(['RMS Error: ', num2str(RMS_error)]);

Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!