Solving Ax=b with both unknowns in A and b

Anyone please help me with this problem. I only know some basic when it comes to matlab and having lots of difficulties trying to code this matrix problem Ax=b that have unknown in both A and b.

The linear system obtained in the form of

where are constants, m is time step, , (initial and boundary values - W at t=0 and W at the last node and outside =0)

How to create the tridiagonal matrix A and vector b in matalb please? I tried to search similiar codes but only found the equations where A and b are known...in this case A and b have symbolic arguments in power form and they are related to the solution

[Torsten]

What is W_{0}^{m} and W_{N+1}^{m} if A is an (NxN) matrix ?

Since you have N+2 unknows, but only N equations, you will need to prescribe the "boundary values" W_{0}^{m} and W_{N+1}^{m}.

Shouldn't b_j be b_i ?

[Walter Roberson]

You can build A using three diag() added together

[Hai Nguyen]

Please correct me if I'm wrong... A and b are extracted from this derived linear governing equation:

is initial condition where time m=0.

Ah yes, it should be

[Hai Nguyen]

Mr Roberson, could you please elaborate a bit more?

[Walter Roberson]

diag(vector1, -1) + diag(vector2) + diag(vector3, 1)

builds a tridiagonal matrix.

[Walter Roberson]

All of the j=i-1 entries are constructed the same way and they all fall along the diagonal just to the left of the major diagonal

Answer 1

[Torsten]

Now if W_{0}^{m} and W_{N+1}^{m} are 0 (I think you meant W_{0}^{m} instead of W_{N}^{m}), you have a system of N equations in W_{1}^{m},...,W_{N}^{m}.

You can use Walter's suggestion to solve for W_{1}^{m+1},...,W_{N}^{m+1}.

I don't understand why you write that there are unknowns in A and b. The W's in A and b are the W's from the last time step (m) and thus known. The unkowns are W_{1}^{m+1},...,W_{N}^{m+1}.

[Hai Nguyen]

For example, given total time = 10 seconds, we need to calculate W's at m=1, 2, 3...10. We know: (m=0)

[Hai Nguyen]

I'm still very new to Matlab and when looking at some similar codes in this problem of Ax=b, I see that A and b are specified. In this case A and b have unknown W's, the unknown x of the main problem is () (). In order to use the solver x=A\b, I think I need to find a way to create A and b symbolically. Is it right?

Actually in this problem, what I really need to calculate are those W's, not . Assumed that N=6 nodes, I need to calculate all (i=1..6) at each time step (m=1..10)

[Torsten]

function main

  C1 = ...;
  C2 = ...;
  C3 = ...;

  m = 25;
  n = 6;

  Wold = zeros(n,1);
  Wnew = zeros(n,1);
  A = zeros(n);
  b = zeros(n,1);

  WNEW = zeros(n,m+1)
  WNEW(:,1) = Wold;

  for k = 1:m
    A(1,1) = C3/C1-4*Wold(1)^3;
    A(1,2) = 2*Wold(2)^3;
    for i = 2:n-1
      A(i,i-1) = 2*Wold(i-1)^3;
      A(i,i) = C3/C1-4*Wold(i)^3;
      A(i,i+1) = 2*Wold(i+1)^3;
    end
    A(n,n-1) = 2*Wold(n-1)^3;
    A(n,n) = C3/C1-4*Wold(n)^3;

    b(1) = -C2/C1-Wold(2)^4+2*Wold(1)^4;
    for i = 2:n-1
      b(i) = -C2/C1-Wold(i+1)^4+2*Wold(i)^4-Wold(i-1)^4;
    end
    b(n) = -C2/C1+2*Wold(n)^4-Wold(n-1)^4;

    dWnew = A\b;
    Wnew = dWnew + Wold;
    Wold = Wnew;

    WNEW(:,k+1) = Wnew;
  end

end

[Hai Nguyen]

Thank you for your great help with the codes. This task might take me a whole year to write.

I ran this code and the result showed W's values are negative and the last node at each column of each W aren't zeros. I think they should be zeros as boundary condition? Could you please explain this to me?

[Torsten]

I specified W_{0}^{m} = W_{7}^{m} = 0 and only saved W_{1}^{m},...,W_{6}^{m} in WNEW.

If you want to set W_{1}^{m} = W_{6}^{m} = 0 and see this in the results, use the code below.

Note that specifying only W_{N}^{m} = 0 doesn't suffice ; you also need a boundary condition on the left (where I took W_{1}^{m} = 0).

If negative values of W don't make sense, you should check again the entries in A and b and the definition of deltaW_{i}^{m}.

function main

  C1 = ...;
  C2 = ...;
  C3 = ...;

  m = 25;
  n = 6;

  Wold = zeros(n,1);
  Wnew = zeros(n,1);
  A = zeros(n);
  b = zeros(n,1);

  WNEW = zeros(n,m+1)
  WNEW(:,1) = Wold;

  for k = 1:m
    A(1,1) = 1.0;
    for i=2:n-1
      A(i,i-1) = 2*Wold(i-1)^3;
      A(i,i) = C3/C1-4*Wold(i)^3;
      A(i,i+1) = 2*Wold(i+1)^3;
    end
    A(n,n) = 1.0;

    b(1) = 0.0;
    for i = 2:n-1
      b(i) = -C2/C1-Wold(i+1)^4+2*Wold(i)^4-Wold(i-1)^4;
    end
    b(n) = 0.0;        

    dWnew = A\b;
    Wnew = dWnew + Wold;
    Wold = Wnew;

    WNEW(:,k+1) = Wnew;
  end

end

[Hai Nguyen]

On the left can you please use the initial condition instead: W_{i}^{0}=0? the problem with one boundary condition and one initial condition also doesn't suffice? In this problem, W_{1} at any time different from zero is to be determined.

The definition of deltaW_{i}^{m}=W_{i}^{m+1} - W_{i}^{m} as edited before.

As this is explicit solving, I guess that it could be the big time step that causes the negative values of W. Can you help re-edit the time step =0.001 instead of 1?

[Torsten]

For further help I need more information about the background of your problem, especially about the underlying continuous equation from which the matrix equation is derived.

[Hai Nguyen]

This problem came from a paper https://doi.org/10.2118/3009-PA

I would like to re-solve the governing equations (coupled of mass conservation and fluid flow equations into a continuity equation) to model hydraulic fracturing process. The solving method is instructed step by step in Apendix A as briefly described below:

(q: flow rate, A: fracture surface area, W: fracture width, h, C_L, q_i: constants)

A-10 is the linear tridiagonal system of the form that I need help to solve using matlab.

[Torsten]

Can you write the original equation (A-1) in W instead of A and q ?

[Hai Nguyen]

In fact, in equation 18 (below) of the paper, the author solved the system in dimensionless form in terms of W_D even though the instruction in Apendix A is for solving (A-1):

Here tau_D is an inverse function of the the maximum length L_D which I don't know how to solve for it, so I treated it as a number.

[Torsten]

Take (A-5) and (A-7) to calculate W_{0}^{m}.

This is the missing boundary value at x=0.

In the last program I posted, you just have to replace b(1) by this value.

[Hai Nguyen]

Thank you Sir for your kindness and patience. If I would like to reduce the time step from 1 to 0.01, which lines of the above program I need to change?

[Torsten]

From your equations, you have to change C1 and C2 accordingly, don't you ?

I wonder where the (delta_x)^2 in the denominator of C1 has gone ...