Help me,master-hands!

[xuej1112]

Main Question or Discussion Point

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 !!!

 

[AiRAVATA]

$$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).$$

 

[xuej1112]

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.

 

[AiRAVATA]

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

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

 

[xuej1112]

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

 

[AiRAVATA]

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:

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.

help diag

Puff!!! Hope this helps.

 

[xuej1112]

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

 

[AiRAVATA]

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

 

[xuej1112]

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

 

[AiRAVATA]

Why are you going backwards if you can go forward?

 

[xuej1112]

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

 

[AiRAVATA]

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 :).

 

[xuej1112]

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