Finite Difference Matrix Help

33 visualizaciones (últimos 30 días)
Justin Yeung
Justin Yeung el 16 de Feb. de 2020
Editada: Fabio Freschi el 21 de Feb. de 2020
So I have a finite difference problem with beam bending. I am trying define a matrix that follows the 4th order ODE for a Central Difference formula.
n= 10 %number of nodes in the beam
%% Step 2: Define the A Matrix
A = (((2*eye(N) +...
diag(ones(N-1,1),1)))+...
diag(ones(N-1,1),-1));
Now this generates a matrix of size N with '2' along the main diagonal. However, a 4th order ODE is different. So I guess my question is, is there a way to add in values '2' off the main diagonal? I've been trying to mess around with the code above, but it keeps saying the "matrix size dimensions must agree.
  2 comentarios
Sindar
Sindar el 16 de Feb. de 2020
check out spdiags. The first example creates the second order second derivative
Srivardhan Gadila
Srivardhan Gadila el 20 de Feb. de 2020
@Justin Yeung are you looking for a matrix which has zeros as diagonal elements, 2's as non-diagonal elements and add this matrix to some other matrix? If not, can you please be more specific or can you give an example?

Iniciar sesión para comentar.

Respuestas (1)

Fabio Freschi
Fabio Freschi el 20 de Feb. de 2020
Editada: Fabio Freschi el 21 de Feb. de 2020
Edit: I changed my answer including a reference and the second order derivative
The coefficients for central differences of different order of accuracy with uniform spacing can be found on wikipedia here.
I assume you need a second order derivative. If it is the case, you can build the matrix using spdiags:
N = 10;
% coefficients (Derivative 2, Accuracy 4 of the wikipedia table)
C = [ones(N,1)/12 4*ones(N,1)/3 -5*ones(N,1)/2 4*ones(N,1)/3 ones(N,1)/12];
% positions along the diagonal
idiag = -2:2;
% matrix
A = spdiags(C,idiag,N,N);
Remember to divide the matrix by the step size dx^2.
The matrix created in this way is sparse (as it is usually done with these problems).

Categorías

Más información sobre Operating on Diagonal Matrices en Help Center y File Exchange.

Productos


Versión

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by