How to solve system of ODEs with unknown initial value?

3 次查看（过去 30 天）
Keidi Zyka
Keidi Zyka 2022-12-20
评论： Keidi Zyka 2022-12-22
Ran in:
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);

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-12-20
编辑：Torsten 2022-12-20
Ran in:
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 个评论
显示 1更早的评论
Torsten
Torsten 2022-12-20
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.

请先登录，再进行评论。

更多回答（1 个）

Bjorn Gustavsson
Bjorn Gustavsson 2022-12-20
编辑：Bjorn Gustavsson 2022-12-20
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 个评论
显示 1更早的评论
Bjorn Gustavsson
Bjorn Gustavsson 2022-12-22
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 2022-12-22
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by