Boundary Value Problem

Question

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)$):

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,
    1}]

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');
    figure(1)
        set(1,'OuterPosition',[0 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
    figure(2)
        set(2,'OuterPosition',[scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
    figure(3)
        set(3,'OuterPosition',[0 0 scrsz(3)/2 scrsz(4)/2])
    figure(4)
        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
            plot(sol.y(4,:),sol.y(3,:))
            axis([-1 1 0 .5])

        figure(2)                         % Fx
            plot(d,sol.parameters(1,:),'o','Markersize',3)
            axis([0 2 -0.2 0.2])
            ylabel('Fx')
            xlabel('d')
            hold on  

        figure(3)                         % Fz
            plot(d,sol.parameters(2,:),'o','Markersize',3)
            axis([0 2 -400 0])
            ylabel('Fz')
            xlabel('d')
            hold on  

        figure(4)                         % M
            plot(d,sol.parameters(3,:),'o','Markersize',3)
            axis([0 2 0 40])
            ylabel('M')
            xlabel('d')
            hold on  

        d = d + solis;
    end
%--------------------------------------------------------------------------
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?

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}]