[Math] Laplace equation in 1D with MATLAB – Dirichlet boundary condition

boundary value problemMATLABordinary differential equations

Here is a Matlab code to solve Laplace 's equation in 1D with Dirichlet's boundary condition u(0)=u(1)=0 using finite difference method

% solve equation -u''(x)=f(x) with the Dirichlet boundary condition 
clear all
close all
N=2;
for j=1:10
a=0;
b=1;

N=2*N;
M(j)=N;
delta_x=(b-a)/N;
k=2;

A=sparse(N-1,N-1);
for i=1:N-1
    if (i==1)
        A(i,i)=2;
        A(i,i+1)=-1;
    elseif(i==N-1)
            A(i,i-1)=-1;
            A(i,i)=2;
    else
        A(i,i-1)=-1;
        A(i,i+1)=-1;
        A(i,i)=2;
    end
end

A=A/((delta_x)^2);


for i=1:N-1
    b(i)=functionf(i*delta_x,k);
end
u=A^(-1)*b';

for i=1:N+1
    u_ex(i)=exact_solution((i-1)*delta_x,k);
end
for i=1:N+1
    if (i==1)
        u_dis(i)=0;
    elseif(i==N+1)
        u_dis(i)=0;
    else
        u_dis(i)=u(i-1,1);
    end
end
for i=1:N+1
    t(i)=(i-1)*delta_x;
end
norm_maxl2(j)=0;
for i=1:N+1
    if (abs(u_dis(i)-u_ex(i)) > norm_maxl2(j))
        norm_maxl2(j)=abs(u_dis(i)-u_ex(i));
    end
end
norm_maxl2(j);
norm_l2(j)=0;
for i=1:N+1
    norm_l2(j)=norm_l2(j)+(u_dis(i)-u_ex(i))^2*delta_x;
end
norm_l2(j)=(norm_l2(j))^(1/2);

norm_maxh1(j)=0;
for i=1:N
    if (abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x) > norm_maxh1(j))
        norm_maxh1(j)=abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x);
    end
end
norm_maxh1(j);

norm_h1(j)=0;
for i=1:N
    norm_h1(j)=norm_h1(j)+(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x)^2*delta_x;
end
norm_h1(j)=(norm_h1(j))^(1/2);

figure
hold on
plot(t,u_ex,'blue', t,u_dis,'red')
ylim([-0.02 0.12]);
hold off
end
%ln(M)
plot(log(M), -log(norm_maxl2),'blue', log(M), -log(norm_l2), 'red', log(M), -log(norm_maxh1), 'cyan', log(M), -log(norm_h1), 'magenta', log(M), 2*log(M),'black')

The file functionf.m

function f=functionf(x,k);

if (k==1)
    f=2;
end

if (k==2)
    f=-6*x+12*x^2;
end

The file exact_solution.m

function u_ex=exact_solution(x,k);
if (k==1)
    u_ex=x*(1-x);
end

if(k==2)
    u_ex=x^3*(1-x);
end

This codes use uniform mesh on [0,1]. But now I have a problem that solve Laplace equation with boundary condition: u(0)=u(1)=0 with singular mesh

enter image description here

How do I modify my code to do this?

Thanks in advanced.

EDITED: Here is my code to solve the problem

% solve equation -u''(x)=f(x) with the Dirichlet boundary condition 
clear all
close all
N=2;
for j=1:10
a=0;
b=1;

N=2*N;
M(j)=N;
k=2;

for i=1:N+1
    x(i)=1-cos((pi/2)*(i-1)/N);
    if i==1
        delta_x(i)=x(i);
    else
    delta_x(i)=x(i)-x(i-1);
    end
end

A=sparse(N,N);
for i=1:N
    if (i==1)
        A(i,i)=2/(delta_x(i))^2;
        A(i,i+1)=-1/(delta_x(i))^2;
    elseif(i==N)
            A(i,i-1)=-1/(delta_x(i))^2;
            A(i,i)=2/(delta_x(i))^2;
    else
        A(i,i-1)=-1/(delta_x(i))^2;
        A(i,i+1)=-1/(delta_x(i))^2;
        A(i,i)=2/(delta_x(i))^2;
    end
end



for i=1:N
    b(i)=functionf(i*delta_x(i),k);
end
u=A^(-1)*b';

for i=1:N+1
    u_ex(i)=exact_solution((i-1)*delta_x(i),k);
end
for i=1:N+1
    if (i==1)
        u_dis(i)=0;
    elseif(i==N+1)
        u_dis(i)=0;
    else
        u_dis(i)=u(i-1,1);
    end
end
for i=1:N+1
    t(i)=(i-1)*delta_x(i);
end
norm_maxl2(j)=0;
for i=1:N+1
    if (abs(u_dis(i)-u_ex(i)) > norm_maxl2(j))
        norm_maxl2(j)=abs(u_dis(i)-u_ex(i));
    end
end
norm_maxl2(j);
norm_l2(j)=0;
for i=1:N+1
    norm_l2(j)=norm_l2(j)+(u_dis(i)-u_ex(i))^2*delta_x(i);
end
norm_l2(j)=(norm_l2(j))^(1/2);

norm_maxh1(j)=0;
for i=1:N
    if (abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i)) > norm_maxh1(j))
        norm_maxh1(j)=abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i));
    end
end
norm_maxh1(j);

norm_h1(j)=0;
for i=1:N
    norm_h1(j)=norm_h1(j)+(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i))^2*delta_x(i);
end
norm_h1(j)=(norm_h1(j))^(1/2);

figure
hold on
plot(t,u_ex,'blue', t,u_dis,'red')
ylim([-0.02 0.12]);
hold off
end
%ln(M)
plot(log(M), -log(norm_maxl2),'blue', log(M), -log(norm_l2), 'red', log(M), -log(norm_maxh1), 'cyan', log(M), -log(norm_h1), 'magenta', log(M), 2*log(M),'black')

But I think it errors because it can't plot u_dis. How must I modify my code?

EDITED (VER $2$): I modifed my code but it look like not ok 🙁 Help me repair it.

% solve equation -u''(x)=f(x) with the Dirichlet boundary condition 
clear all
close all
N=2;
for j=1:10
a=0;
b=1;

N=2*N;
M(j)=N;
k=2;

for i=1:N+1
    x(i)=1-cos((pi/2)*(i-1)/N);
    if i==1
        delta_x(i)=x(i);
    else
    delta_x(i)=x(i)-x(i-1);
    end
end

A=sparse(N,N);
for i=1:N
    if (i==1)
        A(i,i)=2/(delta_x(i+1))^2;
        A(i,i+1)=-1/(delta_x(i+1))^2;
    elseif(i==N)
            A(i,i-1)=-1/(delta_x(i+1))^2;
            A(i,i)=2/(delta_x(i+1))^2;
    else
        A(i,i-1)=-1/(delta_x(i+1))^2;
        A(i,i+1)=-1/(delta_x(i+1))^2;
        A(i,i)=2/(delta_x(i+1))^2;
    end
end



for i=1:N
    b(i)=functionf(x(i),k);
end
u=A^(-1)*b';

for i=1:N+1
    u_ex(i)=exact_solution(x(i),k);
end
for i=1:N+1
    if (i==1)
        u_dis(i)=0;
    elseif(i==N+1)
        u_dis(i)=0;
    else
        u_dis(i)=u(i-1,1);
    end
end
norm_maxl2(j)=0;
for i=1:N+1
    if (abs(u_dis(i)-u_ex(i)) > norm_maxl2(j))
        norm_maxl2(j)=abs(u_dis(i)-u_ex(i));
    end
end
norm_maxl2(j);
norm_l2(j)=0;
for i=1:N
    norm_l2(j)=norm_l2(j)+(u_dis(i)-u_ex(i))^2*delta_x(i+1);
end
norm_l2(j)=(norm_l2(j))^(1/2);

norm_maxh1(j)=0;
for i=1:N-1
    if (abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i+1)) > norm_maxh1(j))
        norm_maxh1(j)=abs(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i+1));
    end
end
norm_maxh1(j);

norm_h1(j)=0;
for i=1:N-1
    norm_h1(j)=norm_h1(j)+(((u_dis(i+1)-u_ex(i+1))-(u_dis(i)-u_ex(i)))/delta_x(i+1))^2*delta_x(i+1);
end
norm_h1(j)=(norm_h1(j))^(1/2);

figure
hold on
plot(x,u_ex,'blue', x,u_dis,'red')
ylim([-0.02 0.12]);
hold off
end
%ln(M)
plot(log(M), -log(norm_maxl2),'blue', log(M), -log(norm_l2), 'red', log(M), -log(norm_maxh1), 'cyan', log(M), -log(norm_h1), 'magenta', log(M), 2*log(M),'black')

Best Answer

Let us discuss here. Firstly, I would like point out that in your formula of $x_i$, the range of subscript should be $i=1\dots N+1$. In addition, RHS is means the right hand side of your numerical scheme.

For your definition of partitions, the partition points is given by the formula of $x_i$. Then the step is given by $h_i=x_{i+1}-x_{i}$ for $i=1\dots N$ (for the sake of simplicity, I use $h_i$ instead of delta_x). By doing this, the matrix A should be reformed in $h_i$, because the original program you wrote is in the case of uniform mesh. What you should do is to adjust the coefficient delta_x to $h_i$. But you must be care of in this case, the step length is depend on step number. What is more, you also should adjust the RHS of your numerical scheme because in your original code, the RHS, say $f(x_i)$, is given in a uniform partition. But in your singular mesh, the RHS should given in the $x_i$ you defined, which is different from uniform mesh.

In summary, just make two vectors for partition points $x_i$ and step length $h_i$. And then modify the code as I told you, I think you can solve the problem.

A sample code:

clear;
N=10;
ax=0;
bx=1;
px=1-cos(pi/2*(0:N)/N);
h=px(2:end)-px(1:end-1);
exau=@(x)x.^2;
RHSf=@(x)-2+zeros(size(x));          %-u''=f
u=zeros(N+1,1);
u(1)=exau(ax);
u(end)=exau(bx);
A=diag(1./h(2:end)+1./h(1:end-1))-diag(1./h(2:end-1),1)-diag(1./h(2:end-1),-1);   
RHS=RHSf(px(2:end-1)').*(h(1:end-1)+h(2:end))'/2;
RHS(1)=RHS(1)+u(1)/h(1);
RHS(end)=RHS(end)+u(end)/h(end);
u(2:end-1)=A\RHS;
plot(px,u-exau(px'));
% hold on;
% plot(linspace(ax,bx,N+1),exau(linspace(ax,bx,N+1)'),'-r');
max(abs(u-exau(px')))

The numerical scheme of 2nd derivative:

Let $x_i$ is the partition point and $h_i=x_{i+1}-x_i$ is the step length. What is more, $u_i$ denote $u(x_i)$ and $u_{i+1/2}$ denote value of $u(x)$ at the middle point of $[x_i,x_{i+1}]$. \begin{equation} u_i''=\frac{u'_{i+1/2}-u'_{i-1/2}}{h_{i+1}/2+h_i/2} \\=\frac{\frac{u_{i+1}-u_i}{h_{i+1}}-\frac{u_i-u_{i-1}}{h_i}}{h_{i+1}/2+h_i/2} \end{equation} Just render it in a neat form, you will get the numerical scheme.

Related Question