I am trying to implement CIRBF on 1D heat eq using crank Nicolson for the j-1 time level. How to plot approximate sol as compared to exact sol.
"clearvars
clc
% parameters
nx=11;
Tf=1;
dt=0.001;
Lx=pi;
dx=Lx/(nx-1);
h=dx;
nt=Tf/dt;
x=(0:dx:Lx);
t=(0:dt:Tf);
u=zeros(nx,nx);
% Initial condition
u(:,1)=sin(2*x)';
%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
x_values=x';
H_bar=[H1,x_values,ones(size(H,1),1)];
H_inv=pinv(H_bar);
% Conversion matrix
C=[H_bar;G_bar];
C_inv=pinv(C);
% 2nd derivative of u at boundaries
f(1,:)=(G1*H_inv)*u(:,1);
f(nx,:)=(G2*H_inv)*u(:,1);
D=[u(:,1);f(1,:);f(nx,:)];
E=C_inv*D;
%weight
W=zeros(1,nx);
for i=1:nx
W(i)=E(i);
end
% find d2u_dx2,u_dx,u
f(2:nx-1,:)=G_a*(C_inv*D);
%Lxx and Bxx
F=G*C_inv;
Bxx=F(:,1:nx);
M=F(:,nx+1);
N=F(:,nx+2);
% Combining M and N into a single matrix B
Lxx=[M,ones(nx,nx-2),N];
% identity & zero matrix
I=eye(nx);
I_a=I*(dt/2);
O=zeros(nx);
% Final Algebric system
A=[Bxx,-Lxx;I,-I_a];
% %rhs
rhs=zeros(nx,nx);
for i=1:nx
for j=2:nx-1
rhs(2:nx-1,1)=u(2:nx-1,1)+(dt/2)*f(2:nx-1,1);
end
end
B=[O;rhs];
%A=[ 1 1 0 3;2 1 -1 1;3 -1 -1 2;-1 2 3 -1];
%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,nx);
y(1,:)=B(1,:)/L(1,1);
for k=2:N
for i=1:nx
y(k,i)=(B(k,i)-L(k,1:k-1)*y(1:k-1,i))/L(k,k);
end
end
% Solve Ux = y for x using backward substitution
x=zeros(N,nx);
x(N,:)=y(N,:)/U(N,N);
for k=N-1:-1:1
for i=1:nx
x(k,i)=(y(k,i)-U(k,k+1:N)*x(k+1:N,i))/U(k,k);
end
end
% Extract the solution for u and f
u_sol=x(1:nx,:);
f_sol=x(nx+1:end,:);
% Create a linearly spaced vector x_range
x_range=linspace(0,Lx,nx);
% Plot u_sol against x_range
figure;
plot(x_range,u_sol);
title('Numerical Solution of u');
xlabel('x');
ylabel('u');
% % Perform LU decomposition on matrix A
% [L,U]=lu(A);
%
% %Solve Ly = B for y using forward substitution
% y=L\B;
%
% % Solve Ux = y for x using backward substitution
% x_sol=U\y;
%
% % Extract the solution for u and f
% u_sol=x_sol(1:nx,:);
% f_sol=x_sol(nx+1:end,:);
%%exact
uexact=zeros(nx,nx);
for j=1:nx
for i=1:nx
% Compute exact solution at (x(i), t(j))
uexact(i,j)=sin(2*x(i))*exp(-4*t(j));
end
end
%plot(x,uexact(:,2))
%Error and RMS
% Error=sum(x(i,j)-uexact(i,j))^2;
% RMS=sqrt(Error/(nx*nx));"
- Badges
- Science topic
- Similar topics
- Filing