Where are the bugs for this ODE finite difference problem that solve using Newton Raphson method?
조회 수: 3 (최근 30일)
이전 댓글 표시
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
채택된 답변
darova
2019년 12월 2일
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
darova
2019년 12월 3일
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)
추가 답변 (1개)
Thiago Henrique Gomes Lobato
2019년 12월 1일
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
참고 항목
카테고리
Help Center 및 File Exchange에서 Symbolic Math Toolbox에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!