# Boundary Value Problem

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?

• 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? Jan 21, 2013 at 14:05
• 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. Jan 21, 2013 at 14:53
• I added it now. Originally I used the picture not the code because all the Subscripts made it hard to read. Jan 21, 2013 at 15:04
• The variables M0,Fi,Fk seem undefined. Apr 18, 2013 at 4:54
• @xslittlegrass, they are to be found. There are 7 unknowns and 7 boundary conditions. Apr 22, 2013 at 14:55

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


• "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. Jan 14, 2014 at 16:17
• Sure @Juris, I would be glad to learn something new. To be honest, I'm not very confident in calculus... Jan 14, 2014 at 17:13