Help me,master-hands!

  • Thread starter xuej1112
  • Start date
  • #1
18
0

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 !!!
 
Last edited by a moderator:

Answers and Replies

  • #2
173
0
[tex]B_j^0 = 1 \qquad \hbox{(initial condition)},[/tex]

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

[tex]\frac{B_{n+1}^i-B_n^i}{h}=0 \qquad \hbox{(right boundary).[/tex]
 
  • #3
18
0
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.
 
  • #4
173
0
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 [itex]u^i[/itex] be the vector of [itex]B_j^i[/itex] 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:
  • #5
18
0
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
 
  • #6
173
0
As I told you before, the convenient way to do it is to express the system of equations in matrix form, i.e

[tex] 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,[/tex]

or in matrix form, let [itex]B^i = (B_2^i,B_3^i,...,B_j^i,...,B_n^i)^T[/itex]

[tex]B^0 = (1,1,...)^T, \qquad i = 0[/tex]

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

[tex]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[/tex]

[tex]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,[/tex]

[tex] B_{n+1}^i = B_n^i,[/tex]

where

[tex]\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}.[/tex]

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 [itex]\Lambda[/itex], I recommend you to check the Matlab function diag.
Code:
help diag
Puff!!! Hope this helps.
 
Last edited:
  • #7
18
0
thank u so mush!but I think my difference method is implicit difference that is:
[tex] 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},[/tex]

so:
[tex]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},[/tex]

maybe this is a equation system
 
  • #8
173
0
Why are you evaluating the PDE in [itex]t_{i+1}[/itex]?
 
  • #9
18
0
because I use the implicit difference method.It is steady.So,maybe this is a equation system.
 
  • #10
173
0
Why are you going backwards if you can go forward?
 
  • #11
18
0
Thank you so much.
by the way,could I find the analyze solution of this equation by matlab?
 
  • #12
173
0
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 :).
 
  • #13
18
0
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
 

Related Threads for: Help me,master-hands!

  • Last Post
Replies
2
Views
13K
Replies
3
Views
4K
Replies
4
Views
9K
Replies
2
Views
1K
Replies
20
Views
2K
Replies
5
Views
5K
Top