How to solve system of ODEs with unknown initial value?

3 次查看(过去 30 天)
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
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 个评论
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 个评论
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 CenterFile Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by