- #1

- 1,782

- 30

The equation is:

[tex]

\partial_{t}u=\partial^{2}_{x}u

[/tex]

I use the discretisation:

[tex]

u_{i+1,j}-u_{i,j}=s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})+s(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1})

[/tex]

Where [itex]s=\delta t/2\delta x[/itex]. This can be rearranged into the following:

[tex]

su_{i+1,j+1}-(2s+1)u_{i+1,j}+su_{i+1,j-1}=-su_{i,j+1}+(2s-1)u_{i,j}-su_{i,j-1}

[/tex]

This is a matrix equation which is tridiagonal and can be solved very easily using matlab's many inbuilt functions.

So suppose I wish to use the following boundary conditions [itex]\partial_{x}u(t,0)=\partial_{x}u(t,L)=0[/itex] how would I go about doing it? I've heard about people going from an [itex]n\times n[/itex] matrix to an [itex](n-2)\times (n-2)[/itex], solving that and just adding in the boundary conditions after that? I tried that method but my solution blew up annoyingly.

Any suggestions on what I am doing wrong? My code is here:

Code:

```
%This program is meant to test out the Crank-Nicolson scheme using a simple
%nonlinear diffusion scheme.
n=2000;m=30;
t=linspace(0,20,n);% set time and distance scales
x=linspace(0,1,m);
dx=x(2)-x(1);dt=t(2)-t(1); %Increments
s=dt/(2*dx^2);%Useful for the solution.
u=zeros(n,m); %set up solution matrix
p=s*ones(1,m-1); q=-(1+2*s)*ones(1,m);
Z=diag(p,1)+diag(p,-1)+diag(q);
A=Z(2:m-1,2:m-1);
%Add in intial condition:
u(1,:)=exp(-5*(x-0.5).^2);
v=zeros(m-2,1);
%Solve the system
for i=2:n-1
%Construct the RHS for the solution
for j=2:m-1
v(j-1)=s*u(i-1,j+1)-(2*s+1)*u(i-1,j)-s*u(i-1,j-1);
end
%Solve for the new time step
w=A\v;
u(i,2:m-1)=w;
u(i,1)=u(i,2);
u(i,end)=u(i,end-1);
end
```