You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to solve a system of third order non-linear differential equations [SOLVED]
33 views (last 30 days)
Show older comments
Hi to everyone,
I'm dealing with a system of third order non-linear differential equations that governs the equilibrium of a circular plate under large displacements hypotesis.
The system is composed by the following equations:
So do anyone have some suggestion on how to solve numerically these equations?
I know that Matlab have the ODE tool, with ode45 or ode113, to solve numerically system of differential equations, but i don't know how to fit the solution used for that problem into my actual problem. I have tried to look in similar solutions, but also in that case i can't find anything that fit well in my case.
If you have any suggestion, you are welcomed!
EDIT: It seems that i was not too much clear before, so i have to specify that all the derivatives are about the radius 'r' and they are not derivatives of time 't'. So also the 'r' in the equations is a variable, not a datum. The solution have to be found on the radius of the plate, not in time domaine. So the solution span have to be something like:
r_span = [0:R];
Better if the integration step is smaller like R/100 or R/1000. Thanks!
EDIT2: I also have the problem to set different boundary conditions on different point of the plate, so i didn't have an unique 'initial condition' because i can express bcs on the edge, where u=w=u'=w'=0, and in the center of the plate where w''=w'''=0. How can i handle this situation of mixed points on where to apply bcs?
Answers (1)
Stephan
on 22 Jun 2019
Edited: Stephan
on 22 Jun 2019
% create a 1.order system and a function handle
syms u(t) w(t) r h nue
eq(1) = diff(u,t,2) == -1/r*diff(u,t)+u/r-(1-nue)/2*r*(diff(w,t))^2-diff(w,t)*diff(w,t,2);
eq(2) = diff(w,t,3) == -1/r*diff(w,t,2)+1/r^2*w+12/h^2*diff(w,t)*(diff(u,t)+nue*u/r+1/2*(diff(w,t))^2);
[Eqs,Vars] = odeToVectorField(eq);
fun = matlabFunction(Eqs,'Vars',{'t','Y','r','h','nue'});
% calculate numerical
r = 4;
h = 2.5;
nue = 0.3;
tspan = [0 10];
initial_values = [0 0 1 0 1];
[tsol,ysol] = ode45(@(t,Y)fun(t,Y,r,h,nue),tspan,initial_values);
plot(tsol, ysol)
40 Comments
Alessandro Olivieri
on 23 Jun 2019
Edited: Alessandro Olivieri
on 23 Jun 2019
Thanks i will try that as soon as i can. I take that equations from the book of Timoshenko so i assumed to be correct (i checked all the passages that lead to the system). Only thing that i think you missed is that the derivatives are not of t (which i assume that you consider time), but on the radius r. It will work if i change all the 't' with 'r'?
EDIT: it didn't work. It give me the error: "The system does not seem to contain a first-order differential equation for the field 'D(D(w))(t)'."
Stephan
on 23 Jun 2019
here is a Workaround:
% create a 1.order system and a function handle
clear all
close
clc
syms u(r) w(r) nue h t
eq(1) = diff(u,r,2) == -1/r*diff(u,r)+u/r-(1-nue)/2*r*(diff(w,r))^2-diff(w,r)*diff(w,r,2);
eq(2) = diff(w,r,3) == -1/r*diff(w,r,2)+1/r^2*w+12/h^2*diff(w,r)*(diff(u,r)+nue*u/r+1/2*(diff(w,r))^2);
eq = subs(eq,r,t);
[eq,vars] = odeToVectorField(eq)
fun = matlabFunction(eq,'Vars',{'t','Y','nue','h'});
strfun = func2str(fun);
strfun = replace(strfun,'t','r');
odefun = str2func(strfun)
% calculate numerical
nue = 0.3;
h = 5;
tspan = [4 0];
initial = [-1 0 0 1 0];
[rsol, ysol] = ode45(@(r,Y)odefun(r,Y,nue,h), tspan, initial);
plot(rsol, ysol)
Alessandro Olivieri
on 24 Jun 2019
Edited: Alessandro Olivieri
on 24 Jun 2019
Why did you put as initial condition -1 and 1 in two points? I need to put zero (clamped edge) to every initial condition, it that a problem?
Anyway it give me this error: "Error using sym/matlabFunction>checkVarsSubset (line 236)
The free variable nue must be included in the 'Vars' value."
Torsten
on 24 Jun 2019
Why did you put as initial condition -1 and 1 in two points? I need to put zero (clamped edge) to every initial condition, it that a problem?
No problem. In this case, I can even give the analytical solution for your problem: u=w=0.
Alessandro Olivieri
on 24 Jun 2019
Ok Thanks i immagine that. Do you know how to solve the error that occur?
Torsten
on 24 Jun 2019
r = 4;
h = 2.5;
nue = 0.3;
tspan = [0 10];
initial_values = [0 0 1 0 1];
fun = @(t,y) [y(2);-1/r*y(2)+y(1)/r-(1-nue)/(2*r)*y(4)^2-y(4)*y(5);y(4);y(5);-1/r*y(5)+1/r^2*y(3)+12/h^2*y(4)*(y(2)+nue*y(1)/r+1/2*y(4)^2)];
[tsol,ysol] = ode45(fun,tspan,initial_values);
plot(tsol, ysol)
Alessandro Olivieri
on 24 Jun 2019
That are derivatives of r not of t, as i say before. I have just tried this way but it seems not to work.
Torsten
on 24 Jun 2019
h = 2.5;
nue = 0.3;
rspan = [1e-5 10];
initial_values = [0 0 1 0 1];
fun = @(r,y) [y(2);-1/r*y(2)+y(1)/r-(1-nue)/(2*r)*y(4)^2-y(4)*y(5);y(4);y(5);-1/r*y(5)+1/r^2*y(3)+12/h^2*y(4)*(y(2)+nue*y(1)/r+1/2*y(4)^2)];
[rsol,ysol] = ode45(fun,tspan,initial_values);
plot(rsol, ysol)
Alessandro Olivieri
on 24 Jun 2019
Edited: Alessandro Olivieri
on 24 Jun 2019
If i put all the initial values to 0 (as it should be), i obtain a solution composed by all zeros. I think that something is still wrong...
I was thinking that on a circular plate the "initial condition" are set on the centre of the plate (where r is zero) or they should be set on the external boundary (where r=R)?
Because it change it all, i would need to set some boundary condition on the edge (clamped edge solution, so u=w=u'=w'=0) and some on the internal boundary (w''=w'''=0). How can i handle this situation? It seems also that i need more then 5 initial condition.
Alessandro Olivieri
on 24 Jun 2019
Edited: Alessandro Olivieri
on 24 Jun 2019
Sorry i didn't understand you in the first place ahahah, excuse me.
Ok now i will try with bvp4c.
If you could provide some suggestion on how to properly set this for my case would be great.
Alessandro Olivieri
on 24 Jun 2019
Your code or the last one? Anyway the problem in using the BVP is that i have to put boundary conditions on both the middle of the plate (r=0) and the edge of the plate (r=R).
I really don't know how to fit this problem using only the explanations given in the Mathworks website for BVP4C.
Torsten
on 24 Jun 2019
function res = bcfcn(ya,yb)
res = [ya(1);ya(2);ya(3);ya(4);yb(5)-value];
end
sets u=u'=w=w'=0 at r=0 and w''=value at r=R.
"bvpfcn" for "bvp4c" can be taken equal to "fun" for "ode45".
Alessandro Olivieri
on 24 Jun 2019
How did the code recognize where to put the bcs in that way? Is ya always and automatically setted for 0 and yb for the end?
Alessandro Olivieri
on 27 Jun 2019
Hey here i am again. I could not work this out. It give me errors on the number of output (errors: too many output arguments). Do you know how to fix it?
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
That's it:
clear all
close
clc
syms u(r) w(r) nu h t
eq(1) = diff(u,r,2) == -1/r*diff(u,r)+u/r-(1-nu)/2*r*(diff(w,r))^2-diff(w,r)*diff(w,r,2);
eq(2) = diff(w,r,3) == -1/r*diff(w,r,2)+1/r^2*w+12/h^2*diff(w,r)*(diff(u,r)+nu*u/r+1/2*(diff(w,r))^2);
eq = subs(eq,r,t);
R = 0.5;
nu = 0.28;
h = 0.01;
tspan = [0 R];
[rsol, ysol] = bvp4c(@(r,Y)bcfcn(r,Y,nu,h), tspan);
plot(rsol, ysol)
function res=bcfcn(ya,yb)
res = [ya(1);ya(2);ya(3);ya(4);yb(5)];
end
As @Torsten said i just change the odefun to the bcfn. Don't know if i have to change other parameters.
Torsten
on 28 Jun 2019
Edited: Torsten
on 28 Jun 2019
h = 2.5;
nue = 0.3;
bvpfcn = @(r,y) [y(2);-1/r*y(2)+y(1)/r-(1-nue)/(2*r)*y(4)^2-y(4)*y(5);y(4);y(5);-1/r*y(5)+1/r^2*y(3)+12/h^2*y(4)*(y(2)+nue*y(1)/r+1/2*y(4)^2)];
bcfcn = @(ya,yb)[ya(1);ya(2);ya(3);ya(4);yb(5)];
guess = @(r)[0 ;0; 1; 0; 1];
rmesh = linspace(1e-5,10,100);
solinit = bvpinit(rmesh, guess);
sol = bvp4c(bvpfcn, bcfcn, solinit);
Please check whether
-(1-nue)/(2*r)*y(4)^2
or
-(1-nue)/2*r*y(4)^2
is correct in the equation for dy(2)/dr.
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
The first expression is correct. Anyway it didn't work, it give me all zeros solution again.
If i put the guess as you put it give to me the error that the jacobian is singular.
I think the error is int the "Solint" because the guess have to be put only on centrain node (i guess) and not on all the node. So if i can make some guess on the boundary i could not make the same prediction for all the plate. That's just my idea, but i don't know if it is correct.
Torsten
on 28 Jun 2019
Anyway it didn't work, it give me all zeros solution again.
If you take the boundary conditions to be zero, you will get the zero solution. That's what I wrote for several times already.
Bjorn Gustavsson
on 28 Jun 2019
Have a look at your differential equations, if you plug your boundary conditions on the outer edge and shoot inwards, how will u and w change and will that solution satisfy your inner boundary conditions?
Alessandro Olivieri
on 28 Jun 2019
If you have a clamped beam the boundary condition that you will impose are all zero on the clamped boundary and no condition on the other edge. So in my case i have to put all zero to my BC. That because the plate is clamped so u=w=u'=w'=0 on the boundary. In the middle of the plate i could say that the rotation are zero (due to simmetry) so w''=0. Moreover the total shear force in the middle will be zero, so w'''=0.
Now, would you say that is not correct to put this boundary?
Maybe i have to check my differential equation system, i just start to think that i miss the component of the external load.
Alessandro Olivieri
on 28 Jun 2019
Ok i figured out that i miss one last component in my equations that now should be like this:
bvpfcn = @(r,y) [y(2);-1/r*y(2)+y(1)/r-(1-nu)/(2*r)*y(4)^2-y(4)*y(5);y(4);y(5);-1/r*y(5)+1/r^2*y(3)+12/h^2*y(4)*(y(2)+nu*y(1)/r+1/2*y(4)^2)+1/D*q*r/2];
So the missing part was in the last equation with q*r/(2*D), where q is the external load (N/mm^2 or MPa if you prefere) r is the variable radius and D is the stiffness of the plate (Nmm or Nm).
Now even if i put this it give me the error of a singular Jacobian.
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
It is the same as before, just changed the bvpfcn:
R=0.5;
h = 0.01;
nu = 0.28;
q = 0.01;
D = 0.0118;
bvpfcn = @(r,y) [y(2);-1/r*y(2)+y(1)/r-(1-nu)/(2*r)*y(4)^2-y(4)*y(5);y(4);y(5);-1/r*y(5)+1/r^2*y(3)+12/h^2*y(4)*(y(2)+nu*y(1)/r+1/2*y(4)^2)+(1/D)*q*r/2];
bcfcn = @(ya,yb)[ya(1);ya(2);ya(3);ya(4);yb(5)];
guess = @(r)[0 ;0; 0; 0; 0];
rmesh = linspace(1e-12,R,100);
solinit = bvpinit(rmesh, guess);
sol = bvp4c(bvpfcn, bcfcn, solinit);
Maybe the error is in the bcfcn? I don't know what the code consider with ya(1) , ya(2) etc. What ya(1) stand for in reality? For y(1), i.e. u in this case?
Bjorn Gustavsson
on 28 Jun 2019
Are you sure that you not just have a divide-by-zero problem at r==0? What happens if you try a range of [eps R]?
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
What do you mean for eps? The solution for r->0 is asymptotical i think, but i use a small number (1e-12) not zero for this reason.
Torsten
on 28 Jun 2019
@Allesandro Olivieri:
y(1) = u,
y(2) = u'
y(3) = w
y(4) = w'
y(5) = w''
ya(1) = u at r=rstart
ya(2) = u' at r=rstart
ya(3) = w at r=rstart
ya(4) = w' at r=rstart
ya(5) = w'' at r=rstart
yb(1) = u at r=rend
yb(2) = u' at r=rend
yb(3) = w at r=rend
yb(4) = w' at r=rend
yb(5) = w'' at r=rend
Now you can form boundary conditions from these variables, e.g. if you want to specify that
u(0) + w(0) = 3,
you have to set one of the five expressions in bcfcn to
ya(1) + ya(3) - 3
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
Ok i guessed right. But i could not impose any condition on w'''? With my previous code i miss a way to define w''' or can i define simply ya(6) in the bcfcn to set it equal to zero? Because you give me 5 bc but i think i need six of them to properly solve the problem (the six bcs that i mention before and in the EDIT of my original post).
Moreover how should i put the guess? Why did you define two sets of 1 for that? Is that important?
Now i can solve the problem with all zeros at the guess function but the solution seems to be terribly wrong (w=0).
Torsten
on 28 Jun 2019
But the value of w''' is given by the differential equation
w''' = ...
You can't set it independently.
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
But it has no mathematical nor physical sense. Why can't I? I used to do that all the time in other situation. Let's say in the theory of the elastic line the differential equation to be solved is something like: y''(x)=M(x)/EJ but i is nonsense to think that i could not give the value for y''(x_0)=0 or whatever in a certain point (that happen in every hinge or simply supported edge).
To me it seems pretty logic to give a given value to w'''(r=0). With the differential equation i set only the function that define the behaviour of w'''(r) not the value of w''' in a given point
Torsten
on 28 Jun 2019
Edited: Torsten
on 28 Jun 2019
I don't know about elastic theory.
What I know is that IVP or BVP solvers only accept initial/boundary conditions up to order (n-1) for a differential equation of order n.
I think you'd better consult a physics forum now because the MATLAB code to your problem has been set up successfully, but you seem to be unsure about the physical model.
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
I'm really confident about the physical problem, the only problem is the matlab code. I think the bvp is not the right way to solve it in matlab and i'm not sure if Matlab could solve it.
And anyway are you saying that you never saw a bc put on one ofe the n-th order derivatives in a differential equation?? Ok you maybe don't know the elastic theory but also in the equation of motion you could put x''=0 if the acceleration is zero somewhere. To me is realy non-sense to don't have the possibility to put condition on all n-th order of a differential equation.
Torsten
on 28 Jun 2019
Ok you maybe don't know the elastic theory but also in the equation of motion you could put x''=0 if the acceleration is zero somewhere.
Of course you can prescribe x'', but in the interior of the region of interest, not at the boundary.
Alessandro Olivieri
on 28 Jun 2019
Edited: Alessandro Olivieri
on 28 Jun 2019
Why not? Stay on the equation of motion, when you drop something the acceleration (aka x'') at time t=0 is fixed!
Maybe you can not do that in Matlab with BVP or ODE, but please don't say that you cannot prescribe a derivative of the same order of the equation on a boundary. It is both physically and mathematically non-sense what you are saying.
Bjorn Gustavsson
on 28 Jun 2019
Edited: Bjorn Gustavsson
on 28 Jun 2019
NO!!! If you drop something the acceleration will be g! Instantly, from the equation of motion:
Alessandro Olivieri
on 29 Jun 2019
Edited: Alessandro Olivieri
on 29 Jun 2019
Exactly it is fixed for t=0! Ok maybe the drop is not the best example, but c'mon guys in physics there are plenty of situation where you could put condition on the boundary for the n-th order derivative of an n-th order differential equation! Almost all differential equation allow that. Motion, wave, elasticity, etc. It have no physical sense to be unable to set a bc on the n-th order derivatives! And neither a mathematical sense, since the problem stay well posed!
Are you coder or mathematicians/physicist? That are basic things of the first year of university (at least basic physics, just little more advanced maths).
Bjorn Gustavsson
on 1 Jul 2019
Alessandro, show me one example, laid out properly where this is done, and not due to symmetry (in your circular geometry) or "material change". I cannot recall seeing that in any of the problems I've encountered. For the heat-equation I've seen fixed boundary temperature, fixed gradiend and fixed heat-flux, but never anything about fixed .
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)