solving PDE problem : Linear Advection diffusion equation problem

5 visualizaciones (últimos 30 días)
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.

Respuesta aceptada

John D'Errico
John D'Errico el 8 de Mayo de 2020
Editada: John D'Errico el 8 de Mayo de 2020
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.
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.
  3 comentarios
Zeynep Toprak
Zeynep Toprak el 8 de Mayo de 2020
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
Zeynep Toprak el 10 de Mayo de 2020
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.

Iniciar sesión para comentar.

Más respuestas (1)

Bjorn Gustavsson
Bjorn Gustavsson el 8 de Mayo de 2020
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

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by