Compact approximation stencils based on integrated flat radial basis functions.Compact local integrated RBF stencils.
1 view (last 30 days)
Show older comments
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)]);
1 Comment
Abhinaya Kennedy
on 29 Aug 2024
Hey! Try mentioning the problem, the error in code or the expected output to help the community resolve your question.
Answers (0)
See Also
Categories
Find more on Surface and Mesh Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!