Symbolic integration, numerical integration, singularity
6 ビュー (過去 30 日間)
Pengpeng Xu 2020 年 12 月 7 日
回答済み: Divija Aleti 2021 年 4 月 21 日
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!
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!