%Experiment 1.1 External
acting force 
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min. step size, 1/2^n, n = ');
%Input the nodal point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the number for first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
ue=(sin(1)/(2*sinh(1)))*sinh(xe)-0.5*sin(xe);
%Solve by Finite Difference Method
for k = s:power
n = 2^k;
h = (b-a)/n;
x = [a+h : h : b-h];
%Initialize the column vectors
vf=(h^2)*sin(x);
%Consturct the matrix for finite difference method
A=diag(-(2+h^2)*ones(n-1,1))+diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
%Solve for approximated solutions
uf=A\vf';
uh(k-s+1) = uf(node/h);
err(k-s+1) = ue(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
end
end
%Extrapolation of the asymptotic expansion
if power >2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(2*k)*uh(k,j+1)-uh(k,j))/(2^(2*k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
end
%Experiment 1.2 Forcing
Function 
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min. step size, 1/2^n, n = ');
%Input the point for point source acted
act=input('Please input the force acts on x = ');
%Input the point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the number for first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
for j = 1:(ne-1)
if xe(j) <= act
ue(j)=(sinh(act-1)*sinh(xe(j)))/sinh(1);
else
ue(j)=(sinh(xe(j)-1)*sinh(act))/sinh(1);
end
end
%Solve the problem by Fourier Series
uF=zeros(1,ne-1);
for i = 1:5000
uF=uF + ((-2*sin(i*pi*act))/(1+(i^2)*pi^2))*sin(i*pi*xe);
end
%Solve by Finite Difference Method
for k = s:power
n = 2^k;
h = (b-a)/n;
m = round(act/h);
x = [a+h:h:b-h];
%Construct the matrix for finite difference method
A=diag(-(8+2*(h^2))*ones(n-1,1)) + diag((4-h^2)*ones(n-2,1),1)
+diag((4-h^2)*ones(n-2,1),-1);
%Initialize the column vectors
vf=zeros(n-1,1);
vf(m)=4*h;
%Solve for solutions
uf=A\vf;
uh(k-s+1)=uf(node/h);
err(k-s+1)=ue(node/he)-uh(k-s+1);
errF(k-s+1)=uF(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
ratioF(k-s)=errF(k-s+1)/errF(k-s);
end
end
%Extrapolation of the asymptotic expansion
if power > 2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(2*k)*uh(k,j+1)-uh(k,j))/(2^(2*k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
errF(1+k,j)=uF(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
ratioF(1+k,j-1)=errF(1+k,j)/errF(1+k,j-1);
end
end
end
end
%Experiment 1.3 Forcing
Function 
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min step size (1/2^n), n = ');
%Input the point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for the exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
for j = 1:(ne-1)
if j <= (ne/2)
ue(j)=(((2*sinh(xe(j))*cosh(0.5))-(sinh(xe(j))+sinh(xe(j)-1)))/sinh(1))-1;
else
ue(j)=(((2*sinh(xe(j)-1)*cosh(0.5))-(sinh(xe(j))+sinh(xe(j)-1)))/sinh(1))+1;
end
end
%Solve by finite difference method
for k = s:power
n = 2^k;
h = (b-a)/n;
x = [a+h:h:b-h];
%Initialize the column vectors
for j = 1:(n-1)
if j < (n/2)
vf(j)=h^2;
else
vf(j)=-(h^2);
end
end
%Consturct the matrix for finite difference method
A=diag(-(2+h^2)*ones(n-1,1))+diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
%Solve for solution
uf=A\vf';
uh(k-s+1)=uf(node/h);
err(k-s+1)=ue(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
end
end
%Extrapolation of the aysmptotic expansion in three steps
if power > 2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(k)*uh(k,j+1)-uh(k,j))/(2^(k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
end
%Experiment 1.3 Forcing
Function
%The function value at the jumped discontinuous point is the average.
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min step size (1/2^n), n = ');
%Input the point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for the exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
for j = 1:(ne-1)
if j <= (ne/2)
ue(j)=(((2*sinh(xe(j))*cosh(0.5))-(sinh(xe(j))+sinh(xe(j)-1)))/sinh(1))-1;
else
ue(j)=(((2*sinh(xe(j)-1)*cosh(0.5))-(sinh(xe(j))+sinh(xe(j)-1)))/sinh(1))+1;
end
end
%Solve for finite difference method
for k = s:power
n = 2^k;
h = (b-a)/n;
x = [a+h:h:b-h];
%Initialize the column vectors
for j = 1:(n-1)
if j < (n/2)
vf(j)=h^2;
else
vf(j)=-(h^2);
end
end
vf(n/2)=0;
%Consturct the matrix for finite difference method
A=diag(-(2+h^2)*ones(n-1,1))+diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
%Solve for solution
uf=A\vf';
uh(k-s+1)=uf(node/h);
err(k-s+1)=ue(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
end
end
%Extrapolation of the aysmptotic expansion in three steps
if power > 2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(2*k)*uh(k,j+1)-uh(k,j))/(2^(2*k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
end
%Experiment 1.4 Forcing
Function 
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min step size (1/2^n), n = ');
%Input the point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for the exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
for j = 1:(ne-1)
if j <= (ne/2)
ue(j)=(((11/4)*cosh(0.5)+2*sinh(0.5)-3)*sinh(xe(j))/sinh(1))-xe(j);
else
ue(j)=((((11/4)*cosh(0.5)-2*sinh(0.5))*sinh(xe(j)-1)-3*sinh(xe(j)))/sinh(1))
+(xe(j)*xe(j)+2);
end
end
%Solve for finite difference method
for k = s:power
n = 2^k;
h = (b-a)/n;
x = [a+h:h:b-h];
%Initialize the column vectors
for j = 1:(n-1)
if j < (n/2)
vf(j)=x(j)*h^2;
else
vf(j)=-x(j)*x(j)*(h^2);
end
end
vf(n/2)=(1/8)*(h^2);
%Consturct the matrix for finite difference method
A=diag(-(2+h^2)*ones(n-1,1))+diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
%Solve for solution
uf=A\vf';
uh(k-s+1)=uf(node/h);
err(k-s+1)=ue(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
end
end
%Extrapolation of the asymptotic expansion in three steps
if power > 2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(2*k)*uh(k,j+1)-uh(k,j))/(2^(2*k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
end
%Experiment 1.5 Forcing
Function 
%Input the value of n to divide internal [0 ,1]
power=input('Please input the min step size (1/2^n), n = ');
%Input the point for error analysis
node=input('Please input the node for error analysis, x = ');
%Initialize the first step size
for i = 2:power
t=node*(2^i)-fix(node*(2^i))-fix(node);
if t == 0
s = i;
break;
end
end
%Initialize the end-point
a = 0;
b = 1;
%Solve for the exact solutions
ne = 2^8;
he = (b-a)/ne;
xe = [a+he:he:b-he];
for j = 1:(ne-1)
if j <= (ne/2)
ue(j)=((3*cosh(0.5)+2*sinh(0.5)-sin(1) + sin(0.5)*cosh(0.5)
+cos(0.5)*sinh(0.5))*sinh(xe(j))/(2*sinh(1)))-xe(j)-1-(sinh(xe(j)-1)/sinh(1));
else
ue(j)=((sinh(xe(j)-1)*(sin(0.5)*cosh(0.5)-cos(0.5)*sinh(0.5)+3*cosh(0.5)-
2*sinh(0.5)-2)-sinh(xe(j))*sin(1))/(2*sinh(1)))+(sin(xe(j))/2);
end
end
%Solve for finite difference method
for k = s:power
n = 2^k;
h = (b-a)/n;
x = [a+h:h:b-h];
%Initialize the column vectors
for j = 1:(n-1)
if j < (n/2)
vf(j)=(x(j)+1)*h^2;
else
vf(j)=-sin(x(j))*(h^2);
end
end
vf(n/2)=((x(n/2)+1-sin(x(n/2)))/2)*(h^2);
%Consturct the matrix for finite difference method
A=diag(-(2+h^2)*ones(n-1,1))+diag(ones(n-2,1),1)+diag(ones(n-2,1),-1);
%Solve for solution
uf=A\vf';
uh(k-s+1)=uf(node/h);
err(k-s+1)=ue(node/he)-uh(k-s+1);
if k >= s+1
ratio(k-s)=err(k-s+1)/err(k-s);
end
end
%Extrapolation of the asymptotic expansion in three steps
if power > 2
switch power
case 3, loop = 1;
case 4, loop = 2;
otherwise, loop = 3;
end
for k = 1:loop
for j = 1:(power-s+1-k)
uh(1+k,j)=(2^(2*k)*uh(k,j+1)-uh(k,j))/(2^(2*k)-1);
err(1+k,j)=ue(node/he)-uh(1+k,j);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
end
%Experiment 2.1 Forcing
Function 
% Find the Exact Solution By Fourier Series
%Initialize the end-point
a = 0;b = 1;
%Initialize the step sizeh1 = (b-a)/2^3;
%Meshgrid the domain with the given step size
[x,y]=meshgrid(a:h1:b,a:h1:b);ue=zeros(2^3+1);
%Summation for Fourier Series
for i = 1:1000 for j = 1:1000 aij=(-4*sin(i*pi*actx)*sin(j*pi*acty))/((i^2+j^2)*(pi)^2); ue=ue+aij*sin(i*pi*x).*sin(j*pi*y); endend
%Save the solution in a file
save ue1 ue -ascii -double
%Experiment 2.1 Forcing
Function
%Initialize the min. step size
power=input('Please input the min. step size, 1/2^n, n = ');
%Input the nodal point for error analysis
errx=input('Please input the point for error analysis x = ');
erry=input(' y = ');
%Initialize the end point
a = 0;
b = 1;
%Loading the exact solution from file
load ue1;
%Initialize the first step size
for i = 2:power
tx=errx*(2^i)-fix(errx*(2^i))-fix(errx);
ty=erry*(2^i)-fix(erry*(2^i))-fix(erry);
if (tx == 0)&(ty == 0)
s = i;
break;
end
end
%Solve by Finite Difference Method
for m = s:power
n = 2^m;
h = (b-a)/n;
node = (n-1)^2;
upone = ones(node,1);
for k = 1:(n-2)
upone(k*(n-1))=0;
end
A1=spdiags(-4*ones(node,1),0,node,node)+spdiags(ones(node,1),-(n-1),node,node)+
spdiags(ones(node,1),(n-1),node,node);
A4=spdiags(upone,-1,node,node);
A5=rot90(A4,2);
A=A1+A4+A5;
f=zeros(node,1);
f((node+1)/2)=-1;
uh=A\f;
ud=zeros(n+1);
for i = 1:(n-1)
for j = 1:(n-1)
ud(i+1,j+1)=uh((i-1)*(n-1)+j);
end
end
u(m-s+1)=ud((n*erry)+1,(n*errx)+1);
err(m-s+1)=u(m-s+1)-ue1(8*erry+1,8*errx+1);
if m >= s+1
ratio(m-s)=err(m-s+1)/err(m-s);
end
end
%Extrapolation of the asymptotic expansion
if power > 2
switch power
case 3, loop=1;
case 4, loop=2;
otherwise, loop=3;
end
end
for k = 1:loop
for j = 1:(power-s+1-k)
u(1+k,j)=(2^(2*k)*u(k,j+1)-u(k,j))/(2^(2*k)-1);
err(1+k,j)=u(1+k,j)-ue1(8*erry+1,8*errx+1);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
% Experiment 2.2 Forcing
Function 
% Find the Exact Solution By Fourier Series
%Initialize the end-point
a = 0;b = 1;
%Initialize the step sizeh1 = (b-a)/2^3;
%Meshgrid the domain with the given step size
[x,y]=meshgrid(a:h1:b,a:h1:b);ue=zeros(2^3+1);
%Summation for Fourier Series
for i = 1:1000
for j = 1:1000
aij=4*((cos(j*pi)-cos(j*pi/2))*(cos(i*pi/2)-1)+
(cos(j*pi/2)-1)*(cos(i*pi)-cos(i*pi/2)))/(i*j*(i^2+j^2)*(pi)^4); ue=ue+aij*sin(i*pi*x).*sin(j*pi*y); endend
%Save the solution in a file
save ue1 ue -ascii -double
% Experiment 2.2 Forcing
Function
%Function Value at the junction is set to be zero
%Initialize the min. step size
power=input('Please input the min. step size, 1/2^n, n = ');%Input
the nodal point for error analysis
errx=input('Please input the point for error analysis x = ');
erry=input(' y = ');
%Initialize the end pointa = 0;b = 1;%Loading the exact solution from fileload exact2;ue=exact2;
%Initialize the first step size
for i = 2:power tx=errx*(2^i)-fix(errx*(2^i))-fix(errx); ty=erry*(2^i)-fix(erry*(2^i))-fix(erry); if (tx == 0)&(ty == 0) s = i; break; endend
%Solve By Finite Difference method
for m = s:power
n = 2^m;
h = (b-a)/n;
node = (n-1)^2;
upone = ones(node,1);
for k = 1:(n-2)
upone(k*(n-1))=0;
end
A1=spdiags(-4*ones(node,1),0,node,node)+spdiags(ones(node,1),-(n-1),node,node)+
spdiags(ones(node,1),(n-1),node,node);
A4=spdiags(upone,-1,node,node);
A5=rot90(A4,2);
A=A1+A4+A5;
f=zeros(node,1);
for i = 1:n-2
if i < n/2
for j = (n/2)+1:(n-1)
f((i-1)*(n-1)+j)=-h^2;
end
else
for j = 1:(n/2)-1
f(i*(n-1)+j)=-h^2;
end
end
end
uh=A\f;
ud=zeros(n+1);
for i = 1:(n-1)
for j = 1:(n-1)
ud(i+1,j+1)=uh((i-1)*(n-1)+j);
end
end
u(m-s+1)=ud((n*erry)+1,(n*errx)+1);
err(m-s+1)=u(m-s+1)-ue(8*erry+1,8*errx+1);
if m >= s+1
ratio(m-s)=err(m-s+1)/err(m-s);
end
end
%Extrapolation of the aysmptotic expansion in three steps
if power > 2 switch power case 3, loop=1;
case 4, loop=2;
otherwise, loop=3;
end
end
for k = 1:loop
for j = 1:(power-s+1-k)
u(1+k,j)=(2^(k)*u(k,j+1)-u(k,j))/(2^(k)-1);
err(1+k,j)=u(1+k,j)-ue(8*erry+1,8*errx+1);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
% Experiment 2.2 Forcing
Function
%Function value at the junction is set to be average
%Initialize the min. step sizepower=input('Please input the min. step size, 1/2^n, n = ');
%Input the nodal point for error analysis
errx=input('Please input the point for error analysis x = ');
erry=input(' y = ');
%Initialize the end point
a = 0;
b = 1;
%Loading the exact solution from file
load exact2;
ue=exact2;
%Initialize the first step size
for i = 2:power
tx=errx*(2^i)-fix(errx*(2^i))-fix(errx);
ty=erry*(2^i)-fix(erry*(2^i))-fix(erry);
if (tx == 0)&(ty == 0)
s = i;
break;
end
end
%Solve By Finite Difference method
for m = s:power
n = 2^m;
h = (b-a)/n;
node = (n-1)^2;
upone = ones(node,1);
for k = 1:(n-2)
upone(k*(n-1))=0;
end
A1=spdiags(-4*ones(node,1),0,node,node)+spdiags(ones(node,1),-(n-1),node,node)+
spdiags(ones(node,1),(n-1),node,node); A4=spdiags(upone,-1,node,node); A5=rot90(A4,2); A=A1+A4+A5; f=zeros(node,1); for i = 1:n-2 if i < n/2 for j = (n/2)+1:(n-1) f((i-1)*(n-1)+j)=-h^2;
end
f((i-1)*(n-1)+n/2)=-0.5*h^2;
else
for j = 1:(n/2)-1
f(i*(n-1)+j)=-h^2;
end
f(i*(n-1)+n/2)=-0.5*h^2;
end
end
for j = 1:(n-1)
f((n/2-1)*(n-1)+j)=-0.5*h^2;
end
uh=A\f;
ud=zeros(n+1);
for i = 1:(n-1)
for j = 1:(n-1)
ud(i+1,j+1)=uh((i-1)*(n-1)+j);
end
end
u(m-s+1)=ud((n*erry)+1,(n*errx)+1);
err(m-s+1)=u(m-s+1)-ue(8*erry+1,8*errx+1);
if m >= s+1
ratio(m-s)=err(m-s+1)/err(m-s);
end
end
%Extrapolation of the aysmptotic expansion in three steps
if power > 2
switch power
case 3, loop=1;
case 4, loop=2;
otherwise, loop=3;
end
end
for k = 1:loop
for j = 1:(power-s+1-k)
u(1+k,j)=(2^(2*k)*u(k,j+1)-u(k,j))/(2^(2*k)-1);
err(1+k,j)=u(1+k,j)-ue(8*erry+1,8*errx+1);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
% Experiment 2.3 Forcing
Function
% Find the Exact Solution By Fourier Series
%Initialize the end-point
a = 0;b = 1;
%Initialize the step sizeh1 = (b-a)/2^3;
%Meshgrid the domain with the given step size
[x,y]=meshgrid(a:h1:b,a:h1:b);
ue=zeros(2^3+1);
%Summation for Fourier Series
for i = 1:1000 for j = 1:1000 p=cos(j*pi)-0.5*cos(j*pi/2)+sin(j*pi/2)/(j*pi);
q=0.5*cos(i*pi/2)-sin(i*pi/2)/(i*pi); r=0.5*cos(j*pi/2)-sin(j*pi/2)/(j*pi);
s=cos(i*pi)-0.25*cos(i*pi/2)+sin(i*pi/2)/(i*pi)-2*(cos(i*pi)-cos(i*pi/2))/(i^2*(pi)^2);
aij=4*(p*q+r*s)/(i*j*(i^2+j^2)*(pi)^4);
ue=ue+aij*sin(i*pi*x).*sin(j*pi*y);
end
end
%Save the solution in a file
save ue1 ue -ascii -double
% Experiment 2.3 Forcing
Function
%Function Value at the junction is set to be zero.
%Initialize the min. step size
power=input('Please input the min. step size, 1/2^n, n = ');%Input the nodal point for error analysis
errx=input('Please input the nodal point for error analysis x = ');
erry=input(' y = ');
%Initialize the end point
a = 0;
b = 1;
%Initialize the first step size
for i = 2:power
tx=errx*(2^i)-fix(errx*(2^i))-fix(errx);
ty=erry*(2^i)-fix(erry*(2^i))-fix(erry);
if (tx == 0)&(ty == 0)
s = i;
break;
end
end
%Loading the exact solution from file
load ue;
%Solve By Finite Difference method
for m = s:power
n = 2^m;
node = (n-1)^2;
h = (b-a)/n;
[x,y]=meshgrid(a:h:b,a:h:b);
upone = ones(node,1);
for k = 1:(n-2)
upone(k*(n-1))=0;
end
A1=spdiags(-4*ones(node,1),0,node,node)+spdiags(ones(node,1),-(n-1),node,node)+
spdiags(ones(node,1),(n-1),node,node);
A4=spdiags(upone,-1,node,node);
A5=rot90(A4,2);
A=A1+A4+A5;
f=zeros(node,1);
for i = 1:n-2
if i < n/2
for j = (n/2)+1:(n-1)
f((i-1)*(n-1)+j)=-x(1,j)^2*y(i+1,j)*h^2;
end
else
for j = 1:(n/2)-1
f(i*(n-1)+j)=-x(1,j+1)*y(i+1,j+1)*h^2;
end
end
end
uh=A\f;
ud=zeros(n+1);
for i = 1:(n-1)
for j = 1:(n-1)
ud(i+1,j+1)=uh((i-1)*(n-1)+j);
end
end
u(m-s+1)=ud((n*erry)+1,(n*errx)+1);
err(m-s+1)=u(m-s+1)-ue(8*erry+1,8*errx+1);
if m >= s+1
ratio(m-s)=err(m-s+1)/err(m-s);
end
end
%Extrapolation of the aysmptotic expansion in three steps
if power > 2
switch power
case 3, loop=1;
case 4, loop=2;
otherwise, loop=3;
end
end
for k = 1:loop
for j = 1:(power-s+1-k)
u(1+k,j)=(2^k*u(k,j+1)-u(k,j))/(2^k-1);
err(1+k,j)=u(1+k,j)-ue(8*erry+1,8*errx+1);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end
% Experiment 2.3 Forcing
Function
%Function Value at the junction is set to be average
%Initialize the min. step size
power=input('Please input the min. step size, 1/2^n, n = ');
%Input the nodal point for error analysis
errx=input('Please input the nodal point for error analysis x = ');
erry=input(' y = ');
%Initialize the end point
a = 0;
b = 1;
%Initialize the first step size
for i = 2:power tx=errx*(2^i)-fix(errx*(2^i))-fix(errx); ty=erry*(2^i)-fix(erry*(2^i))-fix(erry); if (tx == 0)&(ty == 0) s = i; break; endend %Loading the exact solution from file
load ue;
%Solve By Finite Difference method
for m = s:power
n = 2^m;
node = (n-1)^2;
h = (b-a)/n;
[x,y]=meshgrid(a:h:b,a:h:b);
upone = ones(node,1);
for k = 1:(n-2)
upone(k*(n-1))=0;
end
A1=spdiags(-4*ones(node,1),0,node,node)+spdiags(ones(node,1),-(n-1),node,node)+
spdiags(ones(node,1),(n-1),node,node); A4=spdiags(upone,-1,node,node); A5=rot90(A4,2); A=A1+A4+A5; f=zeros(node,1); for i = 1:n-2 if i < n/2 for j = (n/2)+1:(n-1)
f((i-1)*(n-1)+j)=-x(1,j+1)^2*y(i+1,1)*h^2;
end
f((i-1)*(n-1)+n/2)=-0.5*x(1,(n/2)+1)^2*y(i+1,1)*h^2;
else
for j = 1:(n/2)-1
f(i*(n-1)+j)=-x(1,j+1)*y(i+2,1)*h^2;
end
f(i*(n-1)+n/2)=-0.5*x(1,n/2+1)*y(i+2,1)*h^2;
end
end
for j = 1:(n-1)
if j < n/2
f((n/2-1)*(n-1)+j)=-0.5*x(1,j+1)*y(n/2+1,1)*h^2;
elseif j == n/2
f((n/2-1)*(n-1)+j)=-0.5*(x(1,j+1)*y(n/2+1,1)+x(1,j+1)^2*y(n/2+1,1))*h^2;
else
f((n/2-1)*(n-1)+j)=-0.5*x(1,j+1)^2*y(n/2+1,1)*h^2;
end
end
uh=A\f;
ud=zeros(n+1);
for i = 1:(n-1)
for j = 1:(n-1)
ud(i+1,j+1)=uh((i-1)*(n-1)+j);
end
end
u(m-s+1)=ud((n*erry)+1,(n*errx)+1);
err(m-s+1)=u(m-s+1)-ue(8*erry+1,8*errx+1);
if m >= s+1
ratio(m-s)=err(m-s+1)/err(m-s);
end
end
%Extrapolation of the asymptotic expansion in three steps
if power > 2
switch power
case 3, loop=1;
case 4, loop=2;
otherwise, loop=3;
end
end
for k = 1:loop
for j = 1:(power-s+1-k)
u(1+k,j)=(2^(2*k)*u(k,j+1)-u(k,j))/(2^(2*k)-1);
err(1+k,j)=u(1+k,j)-ue(8*erry+1,8*errx+1);
if j >= 2
ratio(1+k,j-1)=err(1+k,j)/err(1+k,j-1);
end
end
end