Seventh order differential equation
4 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Francesco Marchione
il 13 Apr 2023
Modificato: Torsten
il 21 Apr 2023
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):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1355303/image.jpeg)
with these BCs:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1355308/image.jpeg)
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 Commenti
Risposta accettata
Torsten
il 14 Apr 2023
Modificato: Torsten
il 15 Apr 2023
% 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 Commenti
Torsten
il 18 Apr 2023
Modificato: Torsten
il 18 Apr 2023
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])
Più risposte (1)
Torsten
il 13 Apr 2023
Spostato: Torsten
il 13 Apr 2023
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 Commenti
Torsten
il 14 Apr 2023
Look at the examples under
They should show you how to proceed.
If you encounter problems somewhere with your code, you can come back here to ask.
Vedere anche
Categorie
Scopri di più su Ordinary Differential Equations in Help Center e File Exchange
Prodotti
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!