Where are the bugs for this ODE finite difference problem that solve using Newton Raphson method?
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Zhipeng Li
am 1 Dez. 2019
Bearbeitet: Zhipeng Li
am 12 Dez. 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 Kommentare
Akzeptierte Antwort
darova
am 2 Dez. 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 Kommentare
darova
am 3 Dez. 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)
Weitere Antworten (1)
Thiago Henrique Gomes Lobato
am 1 Dez. 2019
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
Siehe auch
Kategorien
Mehr zu Symbolic Math Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!