In solving ordinary differential equations of elasticity, I got the following hint“Unable to find symbolic solution”.How should I modify the solution

syms sigma_theta(r) sigma_r(r) p p0 E mu a b
u=r/E*(-mu*sigma_r+(1-mu)*sigma_theta);
eqn=[2*(sigma_theta-sigma_r)/r==diff(sigma_r,r,1),...
    diff(u,r,1)==1/E*(sigma_r-2*mu*sigma_theta)];
cond=[sigma_theta(a)==-p,sigma_theta(b)==-p0];
S=dsolve(eqn,cond);

Warning: Unable to find symbolic solution.

Answer 1

[Wan Ji]

That means matlab can not find symbolic solution for the ode with the given boundary condintions (Given values at two ends).

Take a look at the general ode, which is always put forward with a given value at only one end called initial conditions.

[Guoyao.Li]

Thanks!

I try to make two equations into one and get the solution as follows.

syms sigma_theta(r) sigma_r(r) p p0 E mu a b
sigma_theta=1/2*r*diff(sigma_r)+sigma_r;
u=r/E*(-mu*sigma_r+(1-mu)*sigma_theta);
eqn=[diff(u,r,1)==1/E*(sigma_r-2*mu*sigma_theta)];
cond=[sigma_r(a)==-p,sigma_r(b)==-p0];
sigma_r=dsolve(eqn,cond);
sigma_theta=1/2*r*diff(sigma_r);
simplify(sigma_theta)

ans = 

simplify(sigma_r)

ans = 

[Wan Ji]

Yes, that worked! This is only because there is a relation between sigma_r and sigma_theta.

I have also found the relation:

syms sigma_theta(r) sigma_r(r) p p0 E mu d1 d2 f(r) c1 c2 s1 s2
syms a b real positive 
% d1 = diff(sigma_r,r,1);
% d2 = diff(sigma_theta,r,1);
u=r/E*(-mu*sigma_r+(1-mu)*sigma_theta);
eqn=[2*(sigma_theta-sigma_r)/r==diff(sigma_r,r,1);...
    diff(u,r,1)==1/E*(sigma_r-2*mu*sigma_theta)]
eq1 = subs(eqn, [diff(sigma_r), diff(sigma_theta)],[d1,d2]);
[d1, d2] = solve(eq1,d1,d2)
d1/d2 % 由于你的sigma_r 和 sigma_theta 对r 对r求导后比例是-2
% 所以令f = sigma_r - sigma_theta
% 有 diff(f,r,1)== -3*f/r 求解就行了
q = dsolve(diff(f,r,1)==(-3*f/r))
% p = int(q/r), get p = c1/(3*r^3)
% p = int(-2*q/r), get p = (2*C1)/(3*r^3)
eq = -s1/(3*r^3) + s2; % 指定两个待定系数；这里的eq就是sigma_theta的表达式
[c1, c2] = solve([subs(eq,r,a)+p; subs(eq, r,b)+p0],[s1,s2]) % 使用边界条件求解两个系数
sigma_theta = subs(eq,[s1,s2],[c1,c2]) % 求得的c1,c2代入eq中
sigma_theta = subs(sigma_theta,[a,b,p0,p],[1,2,5,10]); % 赋予a=1; b=2; p0=5；p=10
rr = 1:0.01:2;
sigma_theta_val = eval(subs(sigma_theta,r,rr));
plot(rr, sigma_theta_val)
xlabel('r'); ylabel('\sigma_\theta')