I have to solve this boundary value problem: $$\frac{\mathrm{d}e_{3x}}{\mathrm{d}l}=(M_0+F_{0z}x-F_{0x}z)e_{3z}$$ $$\frac{\mathrm{d}e_{3z}}{\mathrm{d}l}=-(M_0+F_{0z}x-F_{0x}z)e_{3x}$$ $$\frac{\mathrm{d}x}{\mathrm{d}l}=e_{3x}$$ $$\frac{\mathrm{d}z}{\mathrm{d}l}=e_{3z}$$

With conditions: $$x(0)=x(1)=z(0)=0$$ $$z(1)=1-d$$ $$e_{3x}(0)=0$$ $$e_{3z}(0)=e_{3z}(1)=1$$

$d$ is a parameter in range (0,2), everything else must be computed.

Is this possible in Wolfram Mathematica? The indices confused Mathematica, so I replaced x-related indices with $i$ and z-related indices with $k$. Here is my attempt (I'm trying to find configuration $x(l),z(l)$): enter image description here

And here is the code:

eq1 = D[Subscript[e, i][l], 
   l] == (Subscript[M, 0] + Subscript[F, k] x[l] - 
     Subscript[F, i] z[l]) Subscript[e, k][l]
eq2 = D[Subscript[e, k][l], 
   l] == -(Subscript[M, 0] + Subscript[F, k] x[l] - 
      Subscript[F, i] z[l]) Subscript[e, i][l]
eq3 = D[x[l], l] == Subscript[e, i][l]
eq4 = D[z[l], l] == Subscript[e, k][l]

c1 = x[0] == x[1] == z[0] == 0
c2 = z[1] == 1 - d
c3 = Subscript[e, i][0] == 0
c4 = Subscript[e, k][0] == Subscript[e, k][1] == 1

d = 0.01

s = NDSolve[{eq1, eq2, eq3, eq4, c1, c2, c3, c4}, {x[l], z[l]}, {l, 0,

The physical problem to be solved is this: An elastic rod of length L is placed along $z$ axis. The tips are fixed so they cannot rotate and then the distance betwwen fixators is shortened (they are pushed together) by length D. The rod between becomes curved. Task is to find curvature. $(x(l),z(l))$ are the coords of the rod, $(e_x(l), e_z(l))$ is tangent of the curve and $d=D/L$. $\vec{F}$ and $M$ is the force and torque exerted on the tip of rod upon these deformations.

Here is a MATLAB code that solves this problem of compressing rod and calculates all 7 unknowns for various values of $d$:

function elastigs_stienis

    scrsz = get(0,'ScreenSize');
        set(1,'OuterPosition',[0 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
        set(2,'OuterPosition',[scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
        set(3,'OuterPosition',[0 0 scrsz(3)/2 scrsz(4)/2])
        set(4,'OuterPosition',[scrsz(3)/2 0 scrsz(3)/2 scrsz(4)/2])

    solis = 0.05;  %precizitate
    d = solis;
    parameters = [0 4*pi^2 4*pi*sqrt(d)];  %Fx, Fz, M - minejumi

    sol = bvpinit(linspace(0,1,10),@yinit,parameters,d);

    for i = 1:25
        sol = bvp4c(@odefun,@bcfun,sol,[],d);

        figure(1)                         % stiena konfiguracija
            axis([-1 1 0 .5])

        figure(2)                         % Fx
            axis([0 2 -0.2 0.2])
            hold on  

        figure(3)                         % Fz
            axis([0 2 -400 0])
            hold on  

        figure(4)                         % M
            axis([0 2 0 40])
            hold on  

        d = d + solis;
function y = yinit(l,d)
    y = [0                              %ex
         1                              %ez
         sqrt(d)*(1-cos(2*pi*l))/pi     %x
         l];                            %z
function dydx = odefun(x,y,parameters,d)
    dydx = [(parameters(3)+parameters(2)*y(3)-parameters(1)*y(4))*y(2)  %ex
            -(parameters(3)+parameters(2)*y(3)-parameters(1)*y(4))*y(1) %ez
            y(1)                                                        %x
            y(2)];                                                      %z
function res = bcfun(ya,yb,parameters,d)
    res=[ya(3)           %x(0)=0
         ya(4)           %z(0)=0
         yb(3)           %x(L)=0
         yb(4)-1+d       %z(L)=1-D/L
         ya(1)           %ex(0)=0
         ya(2)-1         %ez(0)=1
         yb(2)-1];       %ez(L)=1

The question remains - can Wolfram Mathematica do the same?

  • $\begingroup$ I was told that in Matlab this requires to state initial configuration and then use solution at $d(x)$ as initial configuration to solve at $d(x+\Delta x)$. Is this also necessary in Mathematica and if so, how to do it? $\endgroup$
    – Džuris
    Jan 21, 2013 at 14:05
  • $\begingroup$ You should provide your code in text form such that those who would want to help you don't have to manually copy down your code. $\endgroup$
    – jVincent
    Jan 21, 2013 at 14:53
  • $\begingroup$ I added it now. Originally I used the picture not the code because all the Subscripts made it hard to read. $\endgroup$
    – Džuris
    Jan 21, 2013 at 15:04
  • $\begingroup$ The variables M0,Fi,Fk seem undefined. $\endgroup$ Apr 18, 2013 at 4:54
  • $\begingroup$ @xslittlegrass, they are to be found. There are 7 unknowns and 7 boundary conditions. $\endgroup$
    – Džuris
    Apr 22, 2013 at 14:55

1 Answer 1


For solving differential equations, you only need as many boundary values as the number of variables of the equation system (multiplied for each variable with the differential order of the variable), 4 in this case. You also need numerical values for your parameters (M0, Fi, Fk) because you solve the equation using NDSolve, which cannot deal with symbolic parameters. Furthermore, you cannot supply more boundary values (initial conditions) than 4 and expect NDSolve to use the excess for solving the parameters - that's just not how differential equations work. Accordingly, here is a working version of your problem:

ClearAll[ei, ek, m0, fi, fk, z, l, x, d];

param = {m0 -> 1., fi -> 1., fk -> 1., d -> 0.01};
eq = {
    ei'[l] == (m0 + fk*x[l] - fi*z[l]) ek[l],
    ek'[l] == -(m0 + fk*x[l] - fi*z[l]) ei[l],
    x'[l] == ei[l],
    z'[l] == ek[l]
init = {
   x[0] == 0,    (* omitting x[1] == 0 *)
   z[1] == 1-d,  (* omitting z[0] == 0 *)
   ei[0] == 0,
   ek[0] == 1    (* omitting ek[1] == 1 *)

sol = NDSolve[Join[eq, init] /. param , {x[l], z[l], ei[l], ek[l]}, {l, 0, 30}];    
Plot[Evaluate[{x[l], z[l], ei[l], ek[l]} /. param /. sol], {l, 0, 30}]

Mathematica graphics

  • $\begingroup$ "that's just not how differential equations work" <-- maybe it's impossible in Mathematica, but the statement of problem is not wrong. I will update my question with description of physical problem and solution in Matlab. $\endgroup$
    – Džuris
    Jan 14, 2014 at 16:17
  • $\begingroup$ Sure @Juris, I would be glad to learn something new. To be honest, I'm not very confident in calculus... $\endgroup$ Jan 14, 2014 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.