Using bvp4c to solve a fourth order nonlinear differential equation with 2 boundary conditions

Hi everybody,

I am struggling with trying to replicate already numerically solved results for a fourth order differential equation (a modified version of the Euler beam equation for two closely spaced free-standing clamped beams):

where E is the youngs modulus (=169GPa), e0 is vacuum permittivity, t is the thickness of each beam (=220nm), V is a voltage, s0 is the initial separation between two closely spaced beams and I:

where w is the width of each beam (230nm).

with boundary conditions :

at x=0 and x=L (where L is the length of the beam - at these points, the beam is fixed).

I've re-arranged the equation so that I have the d4a/dx4 term on the LHS of the equation and C/(s0-2a(x))^2 on the RHS where C=(1/2)*((V^2*e0*t)/E*I) and then tried to solve it numerically using bvp4c as follows:

function mat4bvp
options = bvpset('RelTol', 1e-6, 'AbsTol', 1e-6, 'NMax', 10000)
solinit = bvpinit(linspace(0,5.8e-3,10),[1 0 0 0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit,options);
x = linspace(0,5.8e-3);
y = deval(sol,x);
plot(x,y(1,:));

function dxdy = mat4ode(~,y)
E=169e9; % YMod
w=230e-9; % Width of beam
t=220e-9; % Thickness of beam
I=(t*(w^.3))/12; % Moment of intertia
V=10; % Voltage
e0=8.854187817e-12; % VPerm
const=((0.5*(V.^2)*e0*t)/(E*I)); %All constant values
s0=100e-9; % Initial Separation
dxdy=[y(2)
      y(3)
      y(4)
      const/((s0-2*y(1)).^2)]; %4th order equation
function res = mat4bc(ya,yb)
res=[ya(1) % a=0 at x=0
  ya(2) % da/dx=0 at x=0
  yb(1) % a=0 at x=L
  yb(2)]; % da/dx=0 at x=L

The problem is that I keep getting the error: 'Unable to solve the collocation equations -- a singular Jacobian encountered' or otherwise a mesh limit error and a horrible looking plot. I'm basically aiming to get a positive curve shape that replicates a displacement that starts at 0 for x=0, increases to a max approximately centrally (only to a few nm) and then back to 0 at x=L (which in this case is 5.8um). I seem to have better luck when I try and set all the units so that the y-axis corresponds to nm instead of m, but I'm not actually sure I've set the constants correctly in these units as it becomes slightly difficult when considering e0 and Pa (and V even?). If I want to use m units so that I have a value*e-9 on the y-axis (which means setting the linpsace from 0 to 5.8e-6, and this seems to be where the problems begin..) do I need to increase the tolerance and mesh points for the solver to work? My guesses are somewhat arbitrary as I have no idea what would be a good guess for this system, but does this matter if I have a good enough mesh accuracy and high enough tolerance values?

Any help with this would be greatly appreciated!

~jb311

Answer 1

[Torsten]

My solver solves your ODE without Problems (see the attached file).

Best wishes

Torsten.

[jbradley]

Hi Torsten,

Thank you for taking the time to look over and answer my question. Unfortunately I cannot see any attached file to your post? Apologies if I'm missing something obvious with regards to its location here!

[jbradley]

Have also realised that my linspace should be to e-6 (um) rather than e-3 (mm) in the code I posted above.

Answer 2

[Torsten]

You should at least start with an initial guess that is consistent with the boundary conditions you impose, e.g. [0 0 0 0].

Best wishes

Torsten.

Answer 3

[Torsten]

I calculated the solution for xright = 5.8e-6.

The solution components become so small (in the order of 1e-23 for y(1)) that I guess the Matlab solver identifies the collocation matrix as being singular.

I enclose the results I obtained with my solver for xright = 5.8e-5:

First column: x Second column : y(1) Third column: y(2) Forth column: y(3) Fifth column: y(4)

0.0000000000E+00    0.0000000000E+00    0.0000000000E+00    0.8641164882E-13   -0.8939187932E-08
  
0.1450000000E-05    0.8635498364E-25    0.1160561942E-18    0.7377387736E-13   -0.8492221082E-08
  
0.2900000000E-05    0.3279330972E-24    0.2142574934E-18    0.6178420737E-13   -0.8045254881E-08
  
0.4350000000E-05    0.6995260594E-24    0.2955436441E-18    0.5044263790E-13   -0.7598289327E-08
  
0.5800000000E-05    0.1177288220E-23    0.3608543914E-18    0.3974916801E-13   -0.7151324422E-08
  
0.7250000000E-05    0.1738736559E-23    0.4111294790E-18    0.2970379676E-13   -0.6704360165E-08
  
0.8700000000E-05    0.2362750683E-23    0.4473086492E-18    0.2030652322E-13   -0.6257396555E-08
  
0.1015000000E-04    0.3029572823E-23    0.4703316431E-18    0.1155734643E-13   -0.5810433594E-08
  
0.1160000000E-04    0.3720807837E-23    0.4811382002E-18    0.3456265478E-14   -0.5363471281E-08
  
0.1305000000E-04    0.4419423200E-23    0.4806680589E-18   -0.3996720595E-14   -0.4916509617E-08
  
0.1450000000E-04    0.5109749009E-23    0.4698609559E-18   -0.1080161272E-13   -0.4469548600E-08
  
0.1595000000E-04    0.5777477979E-23    0.4496566269E-18   -0.1695841185E-13   -0.4022588231E-08
  
0.1740000000E-04    0.6409665440E-23    0.4209948060E-18   -0.2246711891E-13   -0.3575628510E-08
  
0.1885000000E-04    0.6994729333E-23    0.3848152261E-18   -0.2732773484E-13   -0.3128669438E-08
  
0.2030000000E-04    0.7522450215E-23    0.3420576186E-18   -0.3154026059E-13   -0.2681711013E-08
  
0.2175000000E-04    0.7983971250E-23    0.2936617136E-18   -0.3510469709E-13   -0.2234753237E-08
  
0.2320000000E-04    0.8371798210E-23    0.2405672397E-18   -0.3802104529E-13   -0.1787796109E-08
  
0.2465000000E-04    0.8679799472E-23    0.1837139244E-18   -0.4028930612E-13   -0.1340839629E-08
  
0.2610000000E-04    0.8903206018E-23    0.1240414937E-18   -0.4190948053E-13   -0.8938837966E-09
  
0.2755000000E-04    0.9038611432E-23    0.6248967223E-19   -0.4288156944E-13   -0.4469286126E-09
  
0.2900000000E-04    0.9083971895E-23   -0.1816794378E-23   -0.4320557382E-13    0.2592334430E-13
  
0.3045000000E-04    0.9038606189E-23   -0.6249325143E-19   -0.4288149458E-13    0.4469798112E-09
  
0.3190000000E-04    0.8903195690E-23   -0.1240449111E-18   -0.4190933268E-13    0.8939330509E-09
  
0.3335000000E-04    0.8679784369E-23   -0.1837170767E-18   -0.4028908905E-13    0.1340885643E-08
  
0.3480000000E-04    0.8371778787E-23   -0.2405700303E-18   -0.3802076463E-13    0.1787837586E-08
  
0.3625000000E-04    0.7983948097E-23   -0.2936640555E-18   -0.3510436036E-13    0.2234788882E-08
  
0.3770000000E-04    0.7522424037E-23   -0.3420594372E-18   -0.3153987718E-13    0.2681739529E-08
  
0.3915000000E-04    0.6994700935E-23   -0.3848164617E-18   -0.2732731604E-13    0.3128689528E-08
  
0.4060000000E-04    0.6409635699E-23   -0.4209954165E-18   -0.2246667786E-13    0.3575638880E-08
  
0.4205000000E-04    0.5777447821E-23   -0.4496565907E-18   -0.1695796360E-13    0.4022587583E-08
  
0.4350000000E-04    0.5109719373E-23   -0.4698602746E-18   -0.1080117419E-13    0.4469535638E-08
  
0.4495000000E-04    0.4419395004E-23   -0.4806667600E-18   -0.3996310563E-14    0.4916483045E-08
  
0.4640000000E-04    0.3720781940E-23   -0.4811363398E-18    0.3456626331E-14    0.5363429804E-08
  
0.4785000000E-04    0.3029549980E-23   -0.4703293086E-18    0.1155763556E-13    0.5810375915E-08
  
0.4930000000E-04    0.2362731497E-23   -0.4473059622E-18    0.2030671617E-13    0.6257321378E-08
  
0.5075000000E-04    0.1738721432E-23   -0.4111265977E-18    0.2970386724E-13    0.6704266193E-08
  
0.5220000000E-04    0.1177277295E-23   -0.3608515136E-18    0.3974908782E-13    0.7151210359E-08
  
0.5365000000E-04    0.6995191621E-24   -0.2955410099E-18    0.5044237697E-13    0.7598153878E-08
  
0.5510000000E-04    0.3279296736E-24   -0.2142553877E-18    0.6178373375E-13    0.8045096748E-08
  
0.5655000000E-04    0.8635403202E-25   -0.1160549498E-18    0.7377315722E-13    0.8492038971E-08
  
0.5800000000E-04   -0.5767269159E-37    0.2744450171E-32    0.8641064645E-13    0.8938980545E-08

Best wishes

Torsten.

[jbradley]

Thanks again for your help Torsten. As matlab has a problem dealing with such small values, would it be better to try and convert the constant values wrt to nm rather than the regular m units? The values for y(1) you have obtained seem much smaller than are quoted in the results I am trying to replicate. Below is the figure I am trying to replicate for differing values of V:

As you can see, y(1) should have ~1nm as it's maximum value for V=10V. Unfortunately, converting everything won't be trivial, but I think it may be the best way round this problem!

~jb311