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

1 回表示 (過去 30 日間)
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 件のコメント
Abhinaya Kennedy
Abhinaya Kennedy 2024 年 8 月 29 日
Hey! Try mentioning the problem, the error in code or the expected output to help the community resolve your question.

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeSurface and Mesh Plots についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by