Where are the bugs for this ODE finite difference problem that solve using Newton Raphson method?
3 views (last 30 days)
Show older comments
Could someone provide me help to solve this Euler-Bernoulli beam equation by using finite difference method and Newton Raphson please.
With boundary value of y(0) = 0 and dy/ds(L) = 0
I am continually getting answers that are nowhere near the results from the bvp4c command.
6 Comments
Accepted Answer
darova
on 2 Dec 2019
Here is my attempt
clc,clear
N = 10; % 10 Set parameters
L = 16*0.3048; %meter
b = 19.625 * 0.0254; d = 1.625 * 0.0254;
p = 26.5851; I = (b*(d^3))/12; h = L/N; F = 0.7436; E = 68.9e6;
alpha = (h^2)/(E*I);
w = p*b*d;
fode = @(s,f) [f(2); -1/E/I*(w*(L-s)+F)*cos(f(1))];
fbound= @(ya,yb) [ya(1)-0; yb(2)-0];
ss = linspace(0,5);
finit = [0 0];
solinit = bvpinit(ss,finit);
sol = bvp4c(fode,fbound,solinit);
plot(sol.x,sol.y(1,:))
7 Comments
darova
on 3 Dec 2019
If f(s) :
then f'(s)
Try the code
syms phi(s) s
syms E I w L F
f = 1/E/I*(w*(L-s)+F)*cos(phi);
df = diff(f,s);
pretty(df)
More Answers (1)
Thiago Henrique Gomes Lobato
on 1 Dec 2019
I didn't check exactly your FEM implementations but one thing that I quickly noticed is that the Newton-Rapson is an iteractive approach, so maybe you get different results simply because your result did not converged. When I iterate about your code I get a very different result that converges actually fast (3 iterations):
% substitute by your images
N = 20 % Set parameters
L = 16*0.3048; %meter
b = 19.625 * 0.0254; d = 1.625 * 0.0254;
p = 26.5851; I = (b*d^3)/12; h = L/N; F = 0.7436; E = 68.9e6;
alpha = h^2/E*I;
w = p*b*d;
S = [h:h:L] ;
y = ones(N,1);
e = ones(N,1);
A = spdiags([e -2*e e],[-1 0 1],N,N);
A(N-1,N) = -2; % fictitious boundaries method
Iterations = 1000;
tol = 1e-20;
for idx =1:Iterations
function_1 = zeros(N,1);
for n = 1:N
function_1(n) = alpha*(w*(L-S(n))+F)*cos(y(n));
end
Fun = A*y + function_1;
% To create the Jacobian of F(y)
Dfdia = zeros(N,1);
for n = 1:N
Dfdia(n) = alpha*(w*(L-S(n))+F)*sin(y(n));
end
J = diag(Dfdia);
step = inv(A+J)*Fun;
y = (y-step);
if norm(step)<tol
Iter = idx
break
end
end
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!