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

4 次查看（过去 30 天）

Zhipeng Li 2019-12-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494122-where-are-the-bugs-for-this-ode-finite-difference-problem-that-solve-using-newton-raphson-method

编辑： Zhipeng Li 2019-12-12

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 个评论
显示 4更早的评论隐藏 4更早的评论

darova 2019-12-2

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

Zhipeng Li 2019-12-2

编辑：Zhipeng Li 2019-12-2

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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

darova 2019-12-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494122-where-are-the-bugs-for-this-ode-finite-difference-problem-that-solve-using-newton-raphson-method#answer_404357

在 MATLAB Online 中打开

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 个评论
显示 5更早的评论隐藏 5更早的评论

darova 2019-12-3

在 MATLAB Online 中打开

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 2019-12-3

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

请先登录，再进行评论。

更多回答（1 个）

Thiago Henrique Gomes Lobato 2019-12-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/494122-where-are-the-bugs-for-this-ode-finite-difference-problem-that-solve-using-newton-raphson-method#answer_404101

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Zhipeng Li 2019-12-2

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.

请先登录，再进行评论。

请先登录，再回答此问题。