How to solve system of ODEs with unknown initial value?

7 visualizzazioni (ultimi 30 giorni)
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);

Risposta accettata

Torsten
Torsten il 20 Dic 2022
Modificato: Torsten il 20 Dic 2022
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
  3 Commenti
Torsten
Torsten il 20 Dic 2022
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
Keidi Zyka il 22 Dic 2022
Torsten thank you very much. It works beautifully.

Accedi per commentare.

Più risposte (1)

Bjorn Gustavsson
Bjorn Gustavsson il 20 Dic 2022
Modificato: Bjorn Gustavsson il 20 Dic 2022
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
  3 Commenti
Bjorn Gustavsson
Bjorn Gustavsson il 22 Dic 2022
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
Keidi Zyka il 22 Dic 2022
@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.

Accedi per commentare.

Tag

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by