Symbolic integration, numerical integration, singularity
4 ビュー (過去 30 日間)
I have a nonlinear differential equation like f = f(x1, x2, x3, a, b, c, r), where x1,x2,x3 are variables that are dependent on time t; a,b,c are parameters which are not going to be set as number but kept as symbols; r is the spatial coordinate, radius actually. What I need to do is int(f*r,r,0,R) to obtain a symbolic differential equation F(x1,x2,x3,a,b,c).
This equation is so long and complex that MATLAB cannot do the symbolic calculation. Numerical integration seems a plausible solution. However, f is singular at r=0 because the expression has 1/r^2..
The code is:
%% problem statement
% W and U are r-dependent, eta and xi are t-dependent, others are parameters
% N can be set as 4, R can be set as 0.05.
lambda0 = zerobess('J', 0 , N);
lambda1 = zerobess('J', 1 , N);
eta = sym('eta',[N,1]);
xi = sym('xi',[N,1]);
W = sym('W',[N,1]); % transverse basis function
for n = 1:N
W(n) = besselj(0,lambda0(n)*r/R);
U = sym('U',[n,1]); % in-plane basis function
for n = 1:N
U(n) = besselj(1,lambda1(n)*r/R); % bessel basis functions
U_0 = T0*(1-nu)/E/delta; % initial strain due to pretension of N0
u_0 = vpa( U_0*r );
w = vpa(dot(eta,W)); % transverse displacement
u = vpa(u_0+dot(xi,U)); % in-plane disp
dw1 = diff(w,r,1);
dw2 = diff(w,r,2);
du1 = diff(u,r,1);
du2 = diff(u,r,2);
eq_w = expand( ...
1 * ( ...
(1-nu^2)/(E*delta) * ( T0*( dw2+1/r*dw1 ) ) ...
+ dw2 * ( du1 + 1/2*dw1^2 + nu*u/r ) ...
+ dw1 * ( du2 + (1+nu)/r*du1 + dw1*dw2 + 1/(2*r)*dw1^2 ) ...
+ rhow * g * (-H+w) ...
eq_u = expand( 1 * ( ...
du2 + 1/r*du1 - u/r^2 + dw1*dw2 + (1-nu)/(2*r)*dw1^2 ...
%% problem happens when perform the Galerkin process
int( eq_w*W*r, r, 0, R)
int( eq_u*U*r, r, 0, R)
回答 (1 件)
Divija Aleti 2021 年 4 月 21 日
The "int" function cannot solve all integrals since symbolic integration is such a complicated task. It is also possible in this case that no analytic or elementary closed-form solution exists.
For more information, have a look at the 'Tips' section of the following link:
Hope this helps!