Seventh order differential equation

Question

Francesco Marchione 2023-4-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1946813-seventh-order-differential-equation

编辑： 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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2023-4-14

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1946813-seventh-order-differential-equation#answer_1215888

编辑：Torsten 2023-4-15

在 MATLAB Online 中打开

% 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更早的评论隐藏 5更早的评论

Torsten 2023-4-17

编辑：Torsten 2023-4-17

在 MATLAB Online 中打开

% Set model parameters

l = 1;

P = 1;

Ga = 1;

Eatilde = 1;

ha = 1;

E1tilde_array = 1:10;

h1 = 1;

E2tilde = 1;

h2 = 1;

xmesh = linspace(-l,l,100);

solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);

hold on

for i = 1:numel(E1tilde_array)

E1tilde = E1tilde_array(i);

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,:);

plot(x,tau)

end

hold off

grid on

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

Torsten 2023-4-18

编辑：Torsten 2023-4-18

在 MATLAB Online 中打开

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 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 ?

请先登录，再进行评论。

Answer 2

Torsten 2023-4-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1946813-seventh-order-differential-equation#answer_1215423

移动：Torsten 2023-4-13

在 MATLAB Online 中打开

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更早的评论隐藏 2更早的评论

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?