Seventh order differential equation

2023 4 13
2 个回答
回答已采纳
更新时间：2023 4 21
7 次查看（30 天）

0 个投票

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:

1 个投票

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

Hello Torsten,
after executing you code here it comes:
Unrecognized function or variable 'bvpfcn'.
Error in @(x,y)bvpfcn(x,y,l,P,Ga,Eatilde,ha,E1tilde,h1,E2tilde,h2)
Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 128)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Torsten
Torsten 2023-4-15
编辑：Torsten 2023-4-15
Did you copy the full code ? It is continued after the plots by the functions "bvpfcn" and "bcfcn".
Hello Torsten, it works!!! Thank you so much.
What if I would like to plot more tau curves in the same graph (varying for example one model parameter such as E1tilde) ?
Torsten
Torsten 2023-4-17
编辑：Torsten 2023-4-17
Ran in:
Ran in:
% 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
Hello,
if I modify the parameters like this:
% Set model parameters
l = 25;
P = 100;
Ga = 1071;
Eatilde = 3000;
ha = 0.3;
E1tilde_array = 1:10;
h1 = 5;
E2tilde = 75000;
h2 = 5;
xmesh = linspace(-l,l,100);
solinit = bvpinit(xmesh, [1 1 1 1 1 1 1 0 0 0]);
the graph is wrong, it should have a parabolic trend
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

0 个投票

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

Torsten
Torsten 2023-4-13
@Walter Roberson comment moved here:
However it is quite valid to set up your questions symbolically, and then to follow the workflow shown in the first example in odeFunction in order to get to a function handle for numeric use.
Torsten
Torsten 2023-4-13
If you have numerical values for all the parameters of your equation, I must correct myself.
You equation is a linear ordinary differential equation of degree 7. Thus numerically solving for the roots of the characteristic polynomial and incorporating the boundary conditions should give you a symbolic solution for it.
So specify the parameters involved, define the equation and boundary conditions and call "dsolve".
@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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by