Where are the bugs for this ODE finite difference problem that solve using Newton Raphson method?

2 views (last 30 days)
Zhipeng Li on 1 Dec 2019
Edited: Zhipeng Li on 12 Dec 2019
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.

Zhipeng Li on 2 Dec 2019
yes, I am trying to find y that denotes as the angle which a tangent to the surface, and s is just denoted as the arc length.
darova on 2 Dec 2019
Shouldn't the curve be horizontal at the end? dy/ds(L) = 0
Zhipeng Li on 2 Dec 2019
Why should it be horizontal, at L is the maxiumum angle, because this is a cantilever beam springboard problem

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,:))

Zhipeng Li on 3 Dec 2019
yes, thats the equation that i am trying to derive from, let it equal to f(y) and try to get df/dy to create a jacobian matrix which y are the angle that i am trying to get. Do I just simply derive cos(y) to (-sin(y)) as I did or i have to derive the 's' that's inside the bracket aswell, like this:
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)
Zhipeng Li on 3 Dec 2019
Yes, finaly got it !! Still not perfect, but I am happy with it.This have been struggled me for days, and again, thank you so much for your help!!

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):
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

1 Comment

Zhipeng Li on 2 Dec 2019
Thank you very much for answering my question. One of the main erro was that I missed a pair of brackts for 'alpha', and now is been edited, but the result is still quite off from the actual solution i got from bvp4c as shown in above comment. I have tried your version of code with bracket bug fixed and is the same, what sort of issue do you think this is? I am not even sure is it math issue or code issue now.