Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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,
    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?

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

1 Answer

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

share|improve this answer
    
"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. –  Juris Jan 14 at 16:17
    
Sure @Juris, I would be glad to learn something new. To be honest, I'm not very confident in calculus... –  István Zachar Jan 14 at 17:13
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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