Hello,
I am having a bit of difficulty with this problem. When I run my code I get an error saying that my vectors are not the same size. Any help is greatly appreciated. The code is provided below.
function ydd=SecDeriv(x,y)
% Calculate the size of the array x
n = length(x);
% Calculate the spacing of the data h
h = x(2)-x(1);
% Second derivative at first point of the array x: 4-point forward difference formula - Error = O(h^2)
ydd(1) = (2*y(1)-5*y(2)+4*y(3)-y(4))/h/h;
% Second derivative at last point of the array x: 4-point backward difference formula - Error = O(h^2)
ydd(n) = (y(n-3)-4*y(n-2)+5*y(n-1)-2*y(n))/h/h;
% Second derivative at the intermediate points: three point central difference formula - Error = O(h^2)
for i=2:n-1
ydd(i) = (y(i-1)-2*y(i)+y(i+1))/h/h;
end
end
This is then used in this script file:
% Plot the exact bending moment and the approximate bending moment
% on the same set of axes
x = [0 24 48 72 96 120 144 168 192 216 240 264 288 312 336 360];
y = [0 -0.111 -0.216 -0.309 -0.386 -0.441 -0.473 -0.479 -0.458 -0.412 -0.345 -0.263 -0.174 -0.090 -0.026 0];
L = 360; E = 29E6; I = 720; q0 = 250;
x1 = linspace(0,L,101);
Mexact = -(q0*L^4/120)*((20/L^2)*(x./L).^3-(12/L^2)*(x./L));
M = E*I*SecDeriv(x,y);
plot(x,M,'*',x1,Mexact,'--');xlabel('x (in)');ylabel('Moment (in-lb)');
legend('Mexp','Mexact');title('Bending Moment');
I am having trouble, I believe, with the line Mexact = -(q0*L^4/120)*((20/L^2)*(x./L).^3-(12/L^2)*(x./L));

 Respuesta aceptada

Star Strider
Star Strider el 26 de Mzo. de 2015

0 votos

Guessing wildly here, but since you want to plot ‘Mexact’ as a function of ‘x1’, you might want to calculate it using ‘x1’ instead of ‘x’:
Mexact = -(q0*L^4/120)*((20/L^2)*(x1./L).^3-(12/L^2)*(x1./L));

2 comentarios

Claire
Claire el 26 de Mzo. de 2015
Yes, you are exactly right. Thank you!
Star Strider
Star Strider el 26 de Mzo. de 2015
My pleasure!

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Christmas / Winter en Centro de ayuda y File Exchange.

Etiquetas

Preguntada:

el 26 de Mzo. de 2015

Comentada:

el 26 de Mzo. de 2015

Community Treasure Hunt

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

Start Hunting!

Translated by