MATLAB Answers

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

11 views (last 30 days)
Zhipeng Li
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.

  6 Comments

Show 3 older comments
Zhipeng Li
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.
Zhipeng Li
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

Sign in to comment.

Accepted Answer

darova
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

Show 4 older comments
Zhipeng Li
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:
Screenshot 2019-12-03 at 20.19.59.png
darova
darova on 3 Dec 2019
If f(s) :
image.png
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
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!!Screenshot 2019-12-03 at 21.49.11.png

Sign in to comment.

More Answers (1)

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

  1 Comment

Zhipeng Li
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by