solving PDE problem : Linear Advection diffusion equation problem

I try to learn how to solve Time dependent PDE in matlab by myself. This question is from Tobin's book.

The initial condtion is

My matlab code is as follows:

n = 100 ; h = 2/n; %n intervals, width 2/n

x = -1 + h* (1: n-1)'; %node locations

D1 = diag(ones(n-1,1),1) - diag(ones(n-1,1),-1) * (1/(2*h)); % the centered finite diff matrix in place of u_x

D2 = toeplitz ([-2 1 zeros(1, n-3)]/ h^2); %diff matrix for u_xx 

f = @(t, u) (0.2*D2*u ) - (D1 * u); % discretized du/dt

u0 = (1 - x.^2) ./ (1 + 50*x.^2); %initial condition

[t, u] = ode15s (f, [0 2], u0); %solve

waterfall (x, t, u)

When I run the code, I get the following error message;

Error using  * 
Inner matrix dimensions must agree.

Error in exercise_713>@(t,u)(0.2*D2*u)-(D1*u)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode15s (line 150)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in exercise_713 (line 16)
[t, u] = ode15s (f, [0 2], u0); %solve

How can I correct this code?

Please give me a feedback whether my matlab code is totally correct or not. Where is the mistake? I try to learn by myself. And I have no one in order to ask. Please give me an idea/comment about my code.

Thanks a lot.

Answer 1

[John D'Errico]

First, it looks as if you are using ode15s for no good reason. Typically, unless you have some initial expectation your IVP is a stiff one, it makes sense to just use ODE45 as the workhorse. ODE45 will be more efficient in general, when the problem is not stiff.

But that is not causing your problem, and you can choose to use any solver you want to use. It should still work.

Your current problem is that you are not thinking about the size of your matrices carefully, and what that implies.

How many unknowns do you have at each time step in t?

Answer: There are n intervals, so n+1 nodes in x. This leaves us with n+1 total unknowns, two of which are implicitly set to zero, the first and last. So in the actual problem you have n-1 unknowns across x to solve for at any time t. Here n = 100, so we expect the length of u to be 99.

What are the sizes of your two matrices D1 and D2? I know this, because you just asked the question to build the matrix D1, and you are using my exact answer here.

https://www.mathworks.com/matlabcentral/answers/523874-write-this-following-matrix?s_tid=prof_contriblnk

Your question there was to build an n by n bi-diagonal matrix D1. It was not to build a matrix of size (n-1) by (n-1). Therefore I told you how to create an nxn matrix.

So your D2 matrix has size 99x99. while D1 has size 100x100.

>> size(D1)
ans =
   100   100
>> size(D2)
ans =
    99    99

An important part of learning to program in MATLAB is to verify simple things about what you are doing. My recommendation is usually to verify every step of your computations. Don't just write a mess of code and then throw it into a solver, hoping the result will make sense.

That comes down to knowing what sizes your matrices will be, and to think about what that implies. A good recommentation of mine is to also test your function that will be passed to the solver. Verify that it runs. Here, that would imply that

f(0,u0)

should produce something besides an error.

If you tried that as a test, to verify the function f was written properly, you would have seen this:

f(0,u0)
Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of
rows in the second matrix. To perform elementwise multiplication, use '.*'.
Error in @(t,u)(0.2*D2*u)-(D1*u) 

And that should have immediately made you check the dimensions of your matrices.

As you are learning, write code SLOWLY. Test EVERYTHING for consistency with your expectations. Eventually you will get to the point that you can write somewhat larger blocks of code without needing to do such basic consistency checks, but you need to learn to walk before you start to run. Even then, the best rule is to ALWAYS go more slowly. Learn to write your code in small blocks. Verify each block. Then when you put it all together you can have some confidence that it MIGHT run the first time.

[Zeynep Toprak]

Thank you for your explanation, and your helping. Then I correct my dimention of the matrix D1 = (diag(ones(n-2,1),1) - diag(ones(n-2,1),-1))/(2*h); by changing with (n-2). Now, both matcies has 99 sizes.

n = 100 ; h = 2/n; %n intervals, width 2/n

x = -1 + h* (1: n-1)'; %node locations

D1 = (diag(ones(n-2,1),1) - diag(ones(n-2,1),-1))/(2*h); % the centered finite diff matrix in place of u_x

D2 = toeplitz ([-2 1 zeros(1, n-3)]/ h^2); %diff matrix for u_xx 

f = @(t, u) (0.2 * D2 *u ) - (D1 * u ); % discretized du/dt

u0 = (1 - x.^2) ./ (1 + 50*x.^2); %initial condition

[t, u] = ode15s (f, [0 2], u0); %solve

waterfall (x, t, u)

I hope, this graph looks correct? What is your opinion about this result?

And I will take a consideration about your suggestion while studying MATLAB.

[Zeynep Toprak]

Also, this question says "solve". But I get only figure. How can I solve it? What is its exact solution? (Note: I see and apply these codes on Tobin's book, thus, I have no further idea on this topic.)

[Zeynep Toprak]

Hello, dear D'Errico, I did very similar question. But I get an error. Why did I get such an error? Can you please explain it? I really want to learn and see my mistake and understand what the error says. Many thanks. My question is here.

Answer 2

[Bjorn Gustavsson]

When you work with problems like this you have to make sure that your matrices end up the sizes you expect them to be. To do that you should make plenty use of size, and whos. When I ran your code I got this:

whos
  Name        Size             Bytes  Class              Attributes

  D1        100x100            80000  double                       
  D2         99x99             78408  double                       
  f           1x1                 32  function_handle              
  h           1x1                  8  double                       
  n           1x1                  8  double                       
  u0         99x1                792  double                       
  x          99x1                792  double                       

So D1*u will not work. To do this you need to make sure that your D1 and D2 matrices are the correct domensions. One way I simplify this is to make the problme small enough for visual inspection - instead of 100 points reduce it to something on the order of 5-10. Then it also becomes more obvious what you'll need to maintain the boundary-conditions too.

This seems to be pretty OK. My 2 additional comments would be to reduce the diffusion-coefficient from 0.2 to 0.05 or something like that, and to compare the solution you get with centred first derivative and an up-stream derivative. Also, as a physicist the boundary conditions make my head hurt, what is supposed to happen at the down-stream boundary?

HTH

[Zeynep Toprak]

and Thank you as well. Your explanation is great! :)