Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help me,master-hands!

  1. May 26, 2009 #1
    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: May 26, 2009
  2. jcsd
  3. May 26, 2009 #2
    [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]
     
  4. May 27, 2009 #3
    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.
     
  5. May 27, 2009 #4
    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 (Text):

    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: May 27, 2009
  6. May 27, 2009 #5
    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
     
  7. May 28, 2009 #6
    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 (Text):

    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 (Text):
    help diag
    Puff!!! Hope this helps.
     
    Last edited: May 28, 2009
  8. May 30, 2009 #7
    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
     
  9. May 31, 2009 #8
    Why are you evaluating the PDE in [itex]t_{i+1}[/itex]?
     
  10. May 31, 2009 #9
    because I use the implicit difference method.It is steady.So,maybe this is a equation system.
     
  11. Jun 1, 2009 #10
    Why are you going backwards if you can go forward?
     
  12. Jun 5, 2009 #11
    Thank you so much.
    by the way,could I find the analyze solution of this equation by matlab?
     
  13. Jun 5, 2009 #12
    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 :).
     
  14. Jul 7, 2009 #13
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Help me,master-hands!
  1. Please help me! (Replies: 20)

  2. Can you help me ? (Replies: 4)

  3. Help me with this (Replies: 2)

Loading...