Seventh order differential equation

10 次查看（过去 30 天）
Francesco Marchione
编辑： Torsten 2023-4-21
Hello,
I would like to solve this system of differential equations in Matlab (and in the end I would like to plot tau and sigma for -l and +l x values):
with these BCs:
where P, h_i, G_i, h_i are numbers (which I would like to define in the code).
Here I started with this:
% y''''''' - a*y'''''' + b*y''' - c*y' = 0
syms s x y(x) Y
Dy = diff(y);
D2y = diff(y,2);
D3y = diff(y,3);
D4y = diff(y,4);
D5y = diff(y,5);
D6y = diff(y,6);
D7y = diff(y,7);
a==10
b==60
c==40
Eqn = D7y - a*D5y + b*D3y -c*Dy == 0;

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-4-14
编辑：Torsten 2023-4-15
Ran in:
% Set model parameters
l = 1;
P = 1;
Ga = 1;
Eatilde = 1;
ha = 1;
E1tilde = 1;
h1 = 1;
E2tilde = 1;
h2 = 1;
xmesh = linspace(-l,l,100);
solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);
sol = bvp4c(@(x,y)bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2), @(ya,yb)bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2),solinit);
x = sol.x;
tau = sol.y(1,:);
sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*sol.y(2,:)- ha/Ga*sol.y(4,:))/(6/(E1tilde*h1^2));
figure(1)
plot(x,tau)
figure(2)
plot(x,sigma)
function dydx = bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*y(2) - ha/Ga*y(4))/(6/(E1tilde*h1^2));
d7ydx7 = Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*y(6) - Eatilde/ha*12/(E1tilde*h1^3)*y(4) + (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*y(2);
dydx = [y(2);y(3);y(4);y(5);y(6);y(7);d7ydx7;y(1);sigma;x*sigma];
end
function res = bcfcn(ya,yb,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
d2sigma_a = ((4/(E1tilde*h1)+2/(E2tilde*h2))*ya(4) - ha/Ga*ya(6))/(6/(E1tilde*h1^2));
d2sigma_b = ((4/(E1tilde*h1)+2/(E2tilde*h2))*yb(4) - ha/Ga*yb(6))/(6/(E1tilde*h1^2));
res = [ya(8);yb(8)+P;ya(9);yb(9);ya(10);yb(10)-P*(h1+ha)/2;d2sigma_a;d2sigma_b;ya(2)-Ga/ha*P/(E1tilde*h1);yb(2)+Ga/ha*2*P/(E2tilde*h2)];
end
  7 个评论
显示 5更早的评论
Torsten
Torsten 2023-4-18
编辑：Torsten 2023-4-18
Try this code whether you get a different result.
It's the analytical solution of your equation.
% Set model parameters
l = 25;
P = 100;
Ga = 1071;
Eatilde = 3000;
ha = 0.3;
E1tilde = 1;
h1 = 5;
E2tilde = 75000;
h2 = 5;
syms x tau(x)
% Solve differential equation
eqn = diff(tau,x,7) - Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*diff(tau,x,5) + Eatilde/ha*12/(E1tilde*h1^3)*diff(tau,x,3) - (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*diff(tau,x) == 0;
tau = dsolve(eqn)
tau0 = tau;
tau1 = diff(tau,x);
tau2 = diff(tau,x,2);
tau3 = diff(tau,x,3);
tau4 = diff(tau,x,4);
tau5 = diff(tau,x,5);
tau6 = diff(tau,x,6);
tau7 = diff(tau,x,7);
sigma = ((4/(E1tilde*h1)+2/(E2tilde*h2))*tau1 - ha/Ga* tau3)/(6/(E1tilde*h1^2));
sigma2 = diff(sigma,x,2);
% Solve for free parameters in solution from boundary conditions
cond1 = int(tau0,x,-1,1) == -P;
cond2 = int(sigma,-l,l) == 0;
cond3 = int(x*sigma,-l,l) == P*(h1+ha)/2;
cond4 = subs(sigma2,x,-l) == 0;
cond5 = subs(sigma2,x,l) == 0;
cond6 = subs(tau1,x,-l) == Ga/ha*P/(E1tilde*h1);
cond7 = subs(tau1,x,l) == -Ga/ha*2*P/(E2tilde*h2);
[A,b] = equationsToMatrix([cond1 cond2 cond3 cond4 cond5 cond6 cond7]);
coeffs = (double(A)\double(b)).';
%Insert boundary conditions in general solution
vars = symvar(tau)
tau0num = subs(tau0,vars(1:7),coeffs);
tau1num = subs(tau1,vars(1:7),coeffs);
tau2num = subs(tau2,vars(1:7),coeffs);
tau3num = subs(tau3,vars(1:7),coeffs);
tau4num = subs(tau4,vars(1:7),coeffs);
tau5num = subs(tau5,vars(1:7),coeffs);
tau6num = subs(tau6,vars(1:7),coeffs);
tau7num = subs(tau7,vars(1:7),coeffs);
sigmanum = subs(sigma,vars(1:7),coeffs);
sigma2num = subs(sigma2,vars(1:7),coeffs);
% Check solution
double(int(tau0num,x,-l,l)+P)
double(int(sigmanum,x,-l,l))
double(int(x*sigmanum,-l,l) - P*(h1+ha)/2)
double(subs(sigma2num,x,-l))
double(subs(sigma2num,x,l))
double(subs(tau1num,x,-l)-Ga/ha*P/(E1tilde*h1))
double(subs(tau1num,x,l)+Ga/ha*2*P/(E2tilde*h2))
error = tau7num - Ga/ha*(4/(E1tilde*h1)+2/(E2tilde*h2))*tau5num + Eatilde/ha*12/(E1tilde*h1^3)*tau3num - (12*Eatilde*Ga/(E1tilde^2*h1^4*ha^2) + 24*Eatilde*Ga/(E1tilde*E2tilde*h1^3*h2*ha^2))*tau1num;
% Plot solution
figure(1)
fplot(error,[-l l])
figure(2)
fplot(tau0num,[-l l])
figure(3)
fplot(sigmanum,[-l l])
Torsten
Torsten 2023-4-21
编辑：Torsten 2023-4-21
@Francesco Marchione
Since I cannot run this code with MATLAB online (it takes too long), I'd be interested whether it gives a different result than the numerical approach. Could you give a short feedback ?

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2023-4-13
移动：Torsten 2023-4-13
A symbolic approach will lead you nowhere because you had to solve for the general roots of a polynomial of degree 7 which is impossible.
So think about a numerical approach.
In order to cope with the integral boundary conditions, I suggest you additionally solve for the functions
F1(y) = integral_{x=-l}^{x=y} tau dx
F2(y) = integral_{x=-l}^{x=y} sigma*x dx
by solving
dF1/dx = tau(x)
dF2/dx = sigma(x)*x
with the boundary conditions
F1(-l) = 0
F1(l) = -P
F2(-l) = 0
F2(l) = P/2 * (h_1+h_a)
Try bvp4c or bvp5c for a solution.
  4 个评论
显示 2更早的评论
Francesco Marchione
Francesco Marchione 2023-4-13
@Torsten yes I can define the values for all the parameters involved. The unknown parameters are just tau and sigma. Can you please help me in writing the code to solve and plot the results?
Torsten
Torsten 2023-4-14
Look at the examples under
https://de.mathworks.com/help/symbolic/dsolve.html
They should show you how to proceed.
If you encounter problems somewhere with your code, you can come back here to ask.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by