*Write Matlab code for solving this BVP using Finite Different Method:
Example of the code:
% Finite difference method
% Approximate the solution of y"=(-2/x)y'+(2/x^2)y+ sin(lnx)/x^2
% for 1<=x<=2 with y(1)=1 and y(2)=2.
p = @(x) (-2/x);
q = @(x) (2/x^2);
r = @(x) (sin(log(x))/x^2);
aa = 1; bb = 2; alpha = 1; beta = 2; n=29; % h = (bb-aa)/(n+1); h=0.1
fprintf(' x w \n');
h = (bb-aa)/(n+1);
a = zeros(1,n+1);
b = zeros(1,n+1);
c = zeros(1,n+1);
d = zeros(1,n+1);
l = zeros(1,n+1);
u = zeros(1,n+1);
z = zeros(1,n+1);
w = zeros(1,n+1);
x = aa+h;
a(1) = 2+h^2*q(x);
b(1) = -1+0.5*h*p(x);
d(1) = -h^2*r(x)+(1+0.5*h*p(x))*alpha;
m = n-1;
for i = 2 : m
x = aa+i*h;
a(i) = 2+h^2*q(x);
b(i) = -1+0.5*h*p(x);
c(i) = -1-0.5*h*p(x);
d(i) = -h^2*r(x);
end
x = bb-h;
a(n) = 2+h^2*q(x);
c(n) = -1-0.5*h*p(x);
d(n) = -h^2*r(x)+(1-0.5*h*p(x))*beta;
l(1) = a(1);
u(1) = b(1)/a(1);
z(1) = d(1)/l(1);
for i = 2 : m
l(i) = a(i)-c(i)*u(i-1);
u(i) = b(i)/l(i);
z(i) = (d(i)-c(i)*z(i-1))/l(i);
end
l(n) = a(n)-c(n)*u(n-1);
z(n) = (d(n)-c(n)*z(n-1))/l(n);
w(n) = z(n);
for j = 1 : m
i = n-j;
w(i) = z(i)-u(i)*w(i+1);
end
i = 0;
fprintf('%5.4f %11.8f\n', aa, alpha);
for i = 1 : n
x = aa+i*h;
fprintf('%5.4f %11.8f\n', x, w(i));
end
i = n+1;
fprintf('%5.4f %11.8f\n', bb, beta);
Matlab code by Thomas algorithm
% To solve the tridiagonal system CX=B,
% where C is a tridiagonal matrix
function X=trisys(A,D,C,B)
%Input: A is the subdiagonal of the coefficient matrix; D is the main diagonal
of the coefficient matrix; C is the superdiagonal of the coefficient matrix; B
is the constant vector of the linear system; Output: X is the solution vector
N=length(B); for k=2:N
mult=A(k-1)/D(k-1);
D(k)=D(k)-mult*C(k-1);
B(k)=B(k)-mult*B(k-1);
end; X(N)=B(N)/D(N);
for k=N-1:-1:1
X(k)=(B(k)-C(k)*X(k+1))/D(k); end
*Write a Matlab program by the method with spatial step size: h=0.1, and time step
size: k=0.002;
Example of the code:
% Program 8.1 Forward difference method for the heat equation with
% D = 1; f(x) = sin2(2_x); u(0; t) = u(1; t) = 0;
% input:space interval [xl,xr], time interval [yb,yt],
% number of space steps M, number of time steps N
% output: solution w
% Example usage: w=heatfd(0,1,0,1,10,250)
%xl=0; xr=1; yb=0, yt=1; M=10; N=250;
function w=heatfd(xl,xr,yb,yt,M,N)
f=@(x) sin(2 _ pi _ x):2;
l=@(t) 0*t; r=@(t) 0*t;
D=1; % diffusion coefficient
h=(xr-xl)/M; k=(yt-yb)/N ;m=M-1; n=N;
sigma=D*k/(h*h); a=diag(1-2*sigma*ones(m,1))+diag(sigma*ones(m-1,1),1);
a=a+diag(sigma*ones(m-1,1),-1);
% (define matrix by built in functions diag() and ones())
lside=l(yb+(0:n)*k);rside=r(yb+(0:n)*k);
w(:,1)=f(xl+(1:m)*h)’; % initial conditions
for j=1:n
w(:,j+1)=a*w(:,j)+sigma*[lside(j);zeros(m-2,1);rside(j)];
end w=[lside;w;rside]; % attach boundary conditions
x=(0:m+1)*h;t=(0:n)*k;
mesh(x,t,w’); % 3-D plot of solution w
%view(60,-120);
view(80,30)
axis([xl xr yb yt -1 1])
Solve the problem by the finite difference method with (space step size) and or small enough
to satisfy the CFL condition; 0.05 h k= h / c
Use Matlab’s mesh command to plot the approximate solutions and the exact solution.
Example of the code:
% Program 8.5 Finite difference solver for 2D Poisson equation
% with Dirichlet boundary conditions on a rectangle
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps
% Output: matrix w holding solution values
% Example usage: w=poisson(0,1,1,2,4,4)
function w=poisson(xl,xr,yb,yt,M,N)
f=@(x,y) 0; % define input function data
g1=@(x) log(x.^2+1); % define boundary values
g2=@(x) log(x.^2+4); % Example 8.8 is shown
g3=@(y) 2*log(y);
g4=@(y) log(y.^2+1);
m=M+1;n=N+1; mn=m*n;
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2;
x=xl+(0:M)*h; % set mesh values
y=yb+(0:N)*k;
A=zeros(mn,mn);b=zeros(mn,1);
for i=2:m-1 % interior points
for j=2:n-1
A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;
A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2;
A(i+(j-1)*m,i+(j-2)*m)=1/k2;A(i+(j-1)*m,i+j*m)=1/k2;
b(i+(j-1)*m)=f(x(i),y(j));
end
end
for i=1:m % bottom and top boundary points
j=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g1(x(i));
j=n;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g2(x(i));
end
for j=2:n-1 % left and right boundary points
i=1;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g3(y(j));
i=m;A(i+(j-1)*m,i+(j-1)*m)=1;b(i+(j-1)*m)=g4(y(j));
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:mn),m,n); %translate from v to w
mesh(x,y,w')
