# Help me,master-hands!

I am working for a week to solve this eqution by matlab,but I am lose.I don't know how to process the boundary condition.

who can tell me some code for solving the equation like this?
View attachment Microsoft Word (2).doc

Thanks a million !!!

Last edited by a moderator:

## Answers and Replies

$$B_j^0 = 1 \qquad \hbox{(initial condition)},$$

$$\frac{B_0^{i+1}-B_0^i}{k}+a\frac{B_1^i - B_0^i}{h} = 0 \qquad \hbox{(left boundary)},$$

$$\frac{B_{n+1}^i-B_n^i}{h}=0 \qquad \hbox{(right boundary).$$

yes, I know these are the boundary conditions.I mean that I don't know how to express them in the matlab program .

but thanks for your reply.

You'll have to specify them in the first time step, then change them according to the evolution of your program. That means modifying them after each step is calculated. Let $u^i$ be the vector of $B_j^i$ entries (j=1,....n), then
Code:
r = 0:h:r*
u = ones(n+1,1) % Initial condition
% Begin time steps
for i=1:m
u(1,1) = u(1,1) - (a*k/h)*(u(2,1)-u(1,1)); % Left boundary condition
u(2:n,1) = Matrix*u(2:n,1); % Matrix form of the PDE
u(n+1,1) = u(n,1); % Right boundary condition
v(:,i) = u(:,1); % This will create a matrix with the value of u^i in each column
end

Last edited:
Thank you so much! you a so kind!

I wrote some code this afternoon.but i think there is some thing wrong

r is the longitudinal coordinates, t is the x-coordinate:

dr = 0.01; %each depth step
Nr = 100; % Choose the number of depth steps
Nt = 40; % Choose the number of time steps
dt = (10*365*24*60*60)/Nt; %Length of each time step in seconds

T = zeros(Nr+1,Nt+1); %Create temperature matrix with Nz+1 rows, and Nt+1 columns

time = [0:120/Nt:120]; %120 months in 10 years
T(:,1) = 1; %Set initial condition

maxiter = 500
for iter = 1:maxiter
%Tlast = T; %Save the last guess
%T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp
for i=2:Nt+1
T(1,i)=(1+dt/dr)*T(1,i-1)-(dt/dr)*T(2,i-1);

for j=2:Nr

depth_2D = (1/(2*dr^2)*dr*(j-1)+(1-dr*(j-1))/(2*dr))*T(j-1,i-1)-(1/(dr^2)*dr*(j-1)+dr*(j-1))*T(j,i-1)+(1/(2*dr^2)*dr*(j-1)+(1-dr*(j-1))/(2*dr))*T(j+1,i-1);

T(j,i) = depth_2D*dt + T(j,i-1);

T(end,i) = T(end-1,i);
end
end
end

As I told you before, the convenient way to do it is to express the system of equations in matrix form, i.e

$$B_j^{i+1} = \frac{\sigma_j k}{h} B_{j-1}^i + \left(1-\frac{2 \sigma_j k}{h^2}-\frac{\mu_j k}{h}-r_j\right)B_j^i + \left(\frac{\sigma_j k}{h^2}+\frac{\mu_j k}{h}\right)B_{j+1}^i,$$

or in matrix form, let $B^i = (B_2^i,B_3^i,...,B_j^i,...,B_n^i)^T$

$$B^0 = (1,1,...)^T, \qquad i = 0$$

$$B_0^{i+1} = B_0^i+\frac{ak}{h} \bigl(B_1^i-B_0^i\bigr),\qquad j=0,\, i=1,...,m$$

$$B_1^{i+1} = \alpha_1 B_0^i + \beta_1 B_1^i + \gamma_1 B_2^i, \qquad j=1,\, i=1,...,m$$

$$B^{i+1} = \left( \begin{array}{cccccc} \alpha_2 & \beta_2 & \gamma_2 & 0 & 0 & \cdots \\ 0 & \alpha_3 & \beta_3 & \gamma_3 & 0 & \cdots \\ 0 & 0 & \alpha_4 & \beta_4 & \gamma_4 & \cdots \\ \vdots & & & & & \end{array}\right) B^i = \Lambda B^i,$$

$$B_{n+1}^i = B_n^i,$$

where

$$\alpha_j = \frac{\sigma_j k}{h},\,\beta_j = 1-\frac{2 \sigma_j k}{h^2}-\frac{\mu_j k}{h}-r_j,\,\gamma_j = \frac{\sigma_j k}{h^2}+\frac{\mu_j k}{h}.$$

Now, the code should look like this:

Code:
u = ones(1,n) ; % Initial condition
MatrixLambda = something % Definition of Lambda.
for i=1:m
u(1,1) = u(1,1) - (a*k/h)*(u(2,1)-u(1,1)); % Left boundary (j = 0)
u(2,1) = a_1 u(1,1) + b_1 u(2,1) + c_1 u(3,1) % j = 1
u(3:n,1) = MatrixLambda*u(3:n,1); % j = 2, ..., n
u(n+1,1) = u(n,1); % Right boundary
V(:,i) = u(1:n+1,1); % Matrix V with the values of u^i in each column
end

To generate the matrix $\Lambda$, I recommend you to check the Matlab function diag.
Code:
help diag

Puff!!! Hope this helps.

Last edited:
thank u so mush!but I think my difference method is implicit difference that is:
$$B_j^i = \frac{\sigma_j k}{h} B_{j-1}^{i+1} + \left(1-\frac{2 \sigma_j k}{h^2}-\frac{\mu_j k}{h}-r_j\right)B_j^{i+1} + \left(\frac{\sigma_j k}{h^2}+\frac{\mu_j k}{h}\right)B_{j+1}^{i+1},$$

so:
$$B^i = \left( \begin{array}{cccccc} \alpha_2 & \beta_2 & \gamma_2 & 0 & 0 & \cdots \\ 0 & \alpha_3 & \beta_3 & \gamma_3 & 0 & \cdots \\ 0 & 0 & \alpha_4 & \beta_4 & \gamma_4 & \cdots \\ \vdots & & & & & \end{array}\right) B^{i+1} = \Lambda B^{i+1},$$

maybe this is a equation system

Why are you evaluating the PDE in $t_{i+1}$?

because I use the implicit difference method.It is steady.So,maybe this is a equation system.

Why are you going backwards if you can go forward?

Thank you so much.
by the way,could I find the analyze solution of this equation by matlab?

What exaclty do you mean?

You want plots and stuff?

I did sort of the same thing for my undergraduate thesis, but my pde was nonlinear, so I had to use the Crank-Nicholson scheme (basically the same thing you did), that's why I was able to help :).

thank you man! I am just back from norway.
I tried the code you give me,but it is not work.
The following is the code,could you help me to check it?

>>dr = 0.005;
Nr = 20;
dt = 0.1;
Nt = 20;
a=0.2339*0.0189;
b=0.2339;
delta=sqrt(0.0073);
u=zeros(Nr+1,Nt+1);
u(1,:)=1;

for j=3:Nr
alpha=(1/2)*delta^2*dr*(j-1);
mu=a-b*dr*(j-1);
rho=dt/(dr)^2;

A=diag((1-2*rho*alpha-(dt/dr)*mu-dt*dr*(j-1))*ones(Nr-2,1))-diag((rho*alpha+(dt/dr)*mu)*ones(Nr-3,1),1)-diag((rho*alpha)*ones(Nr-3,1),-1);
end
MatrixLambda = A; % Definition of Lambda.
for i=1:Nt
u(1,1) = u(1,1) - (a*dt/dr)*(u(2,1)-u(1,1)); % Left boundary (j = 0)
u(2,1) = rho*(1/2)*delta^2*dr*1* u(1,1) + 1-2*rho*(1/2)*delta^2*dr*1-(dt/dr)*(a-b*dr*1)-dt*dr*1* u(2,1) + rho*(1/2)*delta^2*dr*1+(dt/dr)*(a-b*dr*1)* u(3,1); % j = 1
u(3:Nr,1) = MatrixLambda*u(3:Nr,1); % j = 2, ..., n
u(Nr+1,1) = u(Nr,1); % Right boundary
V(:,i) = u(1:Nr+1,1); % Matrix V with the values of u^i in each column
end