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

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.

[darova]

Looks like ode. What about bvp4c?

[Zhipeng Li]

I have found a missing bracket mistake. Now is edited, and the result has been improved, but still not near accurate when compare to bvp4c.

[darova]

Is it graph of ?

[Zhipeng Li]

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]

Shouldn't the curve be horizontal at the end? dy/ds(L) = 0

[Zhipeng Li]

Why should it be horizontal, at L is the maxiumum angle, because this is a cantilever beam springboard problem

Answer 1

[darova]

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]

Thank you so much for your patient and solution, your bvp4c results make my answer looks alot more resonable, since the main curve shape looks similar. I am now suspicious about I have made an error for the Jacobian matrix. how would you differentiate f(y)

[darova]

Plot derivative of f(s)

plot(sol.x,sol.y(1,:))            % f(s)
plot(sol.x,sol.y(2,:))            % f'(s)

[Zhipeng Li]

thank you very much for your help again. your graph shows different result to my jacobian equation, probably that is my main error then.

this would be my last question, dose this diferrential equation looks right to you? or do i have to differentiate the S in the bracket as well?

Dfdia(n) = alpha*(w*(L-S(n))+F)*(-sin(y(n)))

[darova]

This equation

Dfdia(n) = alpha*(w*(L-S(n))+F)*(-sin(y(n)))

is derivative of the original equation?

[Zhipeng Li]

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]

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]

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!!

Answer 2

[Thiago Henrique Gomes Lobato]

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

[Zhipeng Li]

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.