How to solve system of ODEs with unknown initial value?

Hello, I am working to solve a system of ODEs:

,

where

and A being a given 10x10 matrix.

The task requires to calculate the initial value when . So the first step is to find a function .

%Matrix A
A=[zeros(5,5),eye(5,5)
-9	0	4	0	0	-0.09	0	0.04	0	0
0	-4.5	2	0	0	0	-0.045	0.02	0	0
1.333	1.333	-3.666	0	2	0.0133	0.0133	-0.03666	0	0.02
0	0	0	-1.5	0.25	0	0	0	-0.015	0.0025
0	0	1.2	0.2	-3	0	0	0.012	0.002	-0.03];


syms q_01;              %Symbolic scalar.
syms q(t) [10 1];       %Symbolic vector of functions.


 q0=[q_01 4 3 2 1].';   %Initial values for 'displacement'
 dq0=zeros(5,1);        %Initial values for 'speed'
 x0=[q0;dq0];         %Initial values vector

 %System of ODEs
eqns= [diff(q1,t)==A(1,:)*q

    diff(q2,t)==A(2,:)*q 
    diff(q3,t)==A(3,:)*q
    diff(q4,t)==A(4,:)*q
    diff(q5,t)==A(5,:)*q
    diff(q1,t,2)==A(6,:)*q
    diff(q2,t,2)==A(7,:)*q
    diff(q3,t,2)==A(8,:)*q
    diff(q4,t,2)==A(9,:)*q
    diff(q5,t,2)==A(10,:)*q];

 Dq1=diff(q1,t,2);
    Dq2=diff(q2,t,2);
    Dq3=diff(q3,t,2);
    Dq4=diff(q4,t,2);
    Dq5=diff(q5,t,2);

cond =[q1(0)==q_01, q2(0)==4, q3(0)==3, q4(0)==2, q5(0)==1
    Dq1(0)==0, Dq2(0)==0,Dq3(0)==0,Dq4(0)==0,Dq5(0)==0];

 
%using dsolve to receive the solution q1Solt(t) 
[q1Sol(t),q2Sol(t),q3Sol(t),q4Sol(t),q5Sol(t), q6Sol(t),q7Sol(t), q8Sol(t),q59Sol(t), q10Sol(t)]=dsolve(eqns,cond);

Warning: Number of equations greater than number of indeterminates. Trying heuristics to reduce to square system.

Error using mupadengine/feval_internal
Unable to reduce to square system because the number of equations differs from the number of indeterminates.

Error in dsolve>mupadDsolve (line 334)
T = feval_internal(symengine,'symobj::dsolve',sys,x,options);

Error in dsolve (line 203)
sol = mupadDsolve(args, options);

Answer 1

[Torsten]

You have a boundary value problem. I suggest to solve it numerically using BVP4C.

xmesh = linspace(0,3,100);
solinit = bvpinit(xmesh, [1;4;3;2;1;0;0;0;0;0]);
sol = bvp4c(@bvpfcn, @bcfcn, solinit);
format long
sol.y(1,1)

ans = 
  -2.418530283487390

plot(sol.x, sol.y(1,:), '-o')

function dydx = bvpfcn(x,y)
  A=[zeros(5,5),eye(5,5)
     -9	0	4	0	0	-0.09	0	0.04	0	0
     0	-4.5	2	0	0	0	-0.045	0.02	0	0
     1.333	1.333	-3.666	0	2	0.0133	0.0133	-0.03666	0	0.02
     0	0	0	-1.5	0.25	0	0	0	-0.015	0.0025
     0	0	1.2	0.2	-3	0	0	0.012	0.002	-0.03];
  dydx = A*y;
end

function res = bcfcn(ya,yb)
  res = [ya(2)-4;ya(3)-3;ya(4)-2;ya(5)-1;yb(1)-1;ya(6);ya(7);ya(8);ya(9);ya(10)];
end

[Keidi Zyka]

Hey Torsten, thank you for the suggestion. I tried impementing the bvpfc function for this task.

xmesh = linspace(0,3,2^10);
solinit = bvpinit(xmesh, @guess);

sol = bvp4c(@bvpfun, @bcfcn, solinit);

%Coding the equations
function dxdt= bvpfun(t,q)

dxdt=zeros(10,1);

A=[zeros(5,5),eye(5,5)
-9	0	4	0	0	-0.09	0	0.04	0	0
0	-4.5	2	0	0	0	-0.045	0.02	0	0
1.333	1.333	-3.666	0	2	0.0133	0.0133	-0.03666	0	0.02
0	0	0	-1.5	0.25	0	0	0	-0.015	0.0025
0	0	1.2	0.2	-3	0	0	0.012	0.002	-0.03];

dxdt=A*q;
end

The boundary conditions are to be written in the form . I have coded the boundary conditions that include . How can I make use of exactly?

%Setting the boundary conditions
function res= bcfun(ya,yb,q01)
sym q01;
res= [ya(1)-q01
    ya(2)-4
    ya(3)-3
    ya(4)-2
    ya(5)-1
    yb(1)-1
 ];
end

Is function g below a valid function to guess the solution?

function g = guess(x)
g = [sin(x)
     cos(x)];
end

[Torsten]

You need ten boundary conditions for your system. So the condition y1(0)=q01 is superfluous since you have ten conditions: y2,y3,y4,y5,y6,y7,y8,y9 and y10 at x=0 and y1 at x=3. Further, bvp4c doesn't work with symbolic variables.

The guess is not sufficient because you have to prescribe initial values for y1-y10 (thus a 10x1-vector).

Look at the solution I supplied as an answer for more details.

[Keidi Zyka]

Torsten thank you very much. It works beautifully.

Answer 2

[Bjorn Gustavsson]

This works:

% one second-order ODE rephrased as 2 coupled first-order
syms q0 [2 1]
syms a1 a2      
syms q(t) [2 1];
A = [0 1;a1,a2];
Q = dsolve(diff(q,t)==A*q,q(0)==q0);

This one also works:

syms q0 [4 1]   
syms q(t) [4 1];
syms a31 a32 a33 a34 a41 a42 a43 a44
A = [0 0 1 0;0 0 0 1;a31 a32 a33 a34;a41 a42 a43 a44];
%
% A =
% 
% [   0,   0,   1,   0]
% [   0,   0,   0,   1]
% [ a31, a32, a33, a34]
% [ a41, a42, a43, a44] %% Yeah, I didn't put too much work into this one
Q = dsolve(diff(q,t)==A*q,q(0)==q0);

But if you look at the solution to this one you get a couple of screen-fulls of terms for each component. It might be easier if you look at the output from eig (the numerical eigen-values and eigen-vectors), and manually build a solution from the complex exponentials, and then solve for the required initial conditions from that one. (Since A is larger than 4 x 4 we know that "it will be very challenging" to obtain a closed-form analytical roots to the characteristic polynomial).

HTH

[Keidi Zyka]

Thank you for sharing the above script. It motivated to reduce the code to the following

A=[zeros(5,5),eye(5,5)
-9	0	4	0	0	-0.09	0	0.04	0	0
0	-4.5	2	0	0	0	-0.045	0.02	0	0
1.333	1.333	-3.666	0	2	0.0133	0.0133	-0.03666	0	0.02
0	0	0	-1.5	0.25	0	0	0	-0.015	0.0025
0	0	1.2	0.2	-3	0	0	0.012	0.002	-0.03];

syms q_01;
syms q(t) [10 1];

q6(t)=diff(q1,t,2);
q7(t)=diff(q2,t,2);
q8(t)=diff(q3,t,2);
q9(t)=diff(q4,t,2);
q10(t)=diff(q5,t,2);

 q00=[q_01 4 3 2 1].';
 dq00=zeros(5,1);
 x0=[q00;dq00];

 S=dsolve(diff(q,t)==A*q,q(0)==x0);

Nevertheless it is running since a couple of hours...

[Bjorn Gustavsson]

My idea was to use the numerical calculation of the eigen-values and eigen-vectors of your A-matrix. Then use them to build the standard matrix-exponential (provided A is diagonalizable and all that) to get the general solution, and finally use your boundary-conditions to produce a set of algebraic equations for the coefficients. This should work since you have a system of homogenous linear ODEs with constant coefficients. See for example: Matrix_differential_equation.

@Torsten's solution is what most will have to use "for real problems" where the ODE's might have time-variable coefficients, or be nonlinear and no longer homogenous. That I still suggest this approach is because the problem seems to be designed to give you a proper understanding of these ODE-solving steps, and these days I'd be far to lazy to do this and would definitely use the bvp4c-style solution.

[Keidi Zyka]

@Bjorn Gustavsson thank you a lot. The orginal idea of using the nummerical calculation of the EV and EW works just as fine and is a good practice when solving exam exercises by hand.