如何求解初始值未知的ode系统?

10次浏览(过去30天)
Keidi Zyka
Keidi Zyka 2022年12月20日
评论道: Keidi Zyka2022年12月22日
你好,我正在解决一个ode系统:
在哪里
A是一个给定的10x10矩阵。
任务要求计算初始值 .第一步是找到一个函数
%的矩阵
=[0(5、5)、眼睛(5,5)
-9 04 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 02 0.0133 0.0133 -0.03666 0 0.02
00 0 -1.5 0.25 00 0 -0.015 0.0025
00 1.2 0.2 -3 00 0.012 0.002 -0.03];
信谊q_01;%符号标量。
信谊Q (t) [10 1];函数的符号向量。
Q0 =[q_01 4 3 2 1].';位移的初始值
dq0 = 0(5、1);'speed'的初始值
x0 = [q0; dq0];初始值
ode系统
命令= [diff (q1, t) = =(1:) *问
diff (q2, t) = = (2:) * q
diff (q3, t) = = (3:) * q
diff(第四季度,t) = = (4:) * q
diff (q5 t) = = (5:) * q
diff (q1, t, 2) = =(6:) *问
diff (q2, t, 2) = = (7:) * q
diff (q3, t, 2) = = (8:) * q
diff (q4 t 2) = =(9日:)* q
diff (q5 t 2) = = (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);
气孔导度= (q1 (0) = = q_01, q2(0) = = 4,第三季(0)= = 3,第四季度(0)= = 2,q5 (0) = = 1
Dq1 (0) = = 0, Dq2 (0) = = 0, Dq3 (0) = = 0, Dq4 (0) = = 0, Dq5 (0) = = 0);
%使用dsolve接收溶液q1Solt(t)
[q1Sol (t) q2Sol (t) q3Sol (t) q4Sol (t) q5Sol (t) q6Sol (t) q7Sol (t) q8Sol (t) q59Sol (t) q10Sol (t)) = dsolve(方程式,气孔导度);
警告:方程式数目大于未定式数目。尝试启发式化简到平方系统。
使用mupadengine/feval_internal时出错
由于方程的数目不同于不定式的数目,无法化简为平方系统。

解析>mupadDsolve时出现错误(第334行)
T = feval_internal(symengine,'symobj::dsolve',sys,x,options);

解析错误(第203行)
sol = mupadDsolve(参数,选项);

接受的答案

Torsten
Torsten 2022年12月20日
编辑:Torsten 2022年12月20日
你有一个边值问题。我建议用BVP4C数值求解。
Xmesh = linspace(0,3,100);
Solinit = bvpinit(xmesh, [1;4;3;2;1;0;0;0;0]);
Sol = bvp4c(@bvpfcn, @bcfcn, solinit);
格式
sol.y (1, 1)
ans =
-2.418530283487390
情节(溶胶。x, sol.y (1:)“o”
函数Dydx = bvpfcn(x,y)
=[0(5、5)、眼睛(5,5)
-9 04 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 02 0.0133 0.0133 -0.03666 0 0.02
00 0 -1.5 0.25 00 0 -0.015 0.0025
00 1.2 0.2 -3 00 0.012 0.002 -0.03];
dydx = A*y;
结束
函数Res = bcfcn(ya,yb)
res =[丫(2)4,丫(3)3;丫(4)2,丫(5)1;yb(1) 1;丫(6);丫(7);丫(8);丫(9);丫(10)];
结束
3评论
Keidi Zyka
Keidi Zyka 2022年12月22日
Torsten,非常感谢。它工作得很好。

登录评论。

更多答案(1)

Bjorn Gustavsson
Bjorn Gustavsson 2022年12月20日
编辑:Bjorn Gustavsson 2022年12月20日
如此:
% 1次二阶ODE改写为2次耦合一阶ODE
信谊Q0 [2 1]
信谊a1 a2
信谊Q (t) [2 1];
A = [0 1;a1,a2];
Q = dsolve(diff(Q,t)==A* Q, Q (0)==q0);
这个也可以:
信谊Q0 [4 1]
信谊Q (t) [4 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];
% a =
% [0,0,1,0]
% [0,0,0,1]
% [a31, a32, a33, a34]
是啊,我没在这个上面花太多功夫
Q = dsolve(diff(Q,t)==A* Q, Q (0)==q0);
但是,如果你看一下这个问题的解决方案,你会发现每个组件都有满屏的术语。如果查看的输出可能会更简单 eig (数值特征值和特征向量),并手动从复指数建立一个解,然后从那个解求解所需的初始条件。(因为 一个 大于4 × 4我们知道它会是多少 非常 挑战“获得特征多项式的闭形式解析根”)。
HTH
3评论
Keidi Zyka
Keidi Zyka 2022年12月22日
@Bjorn Gustavsson 非常感谢。使用EV和EW的数值计算的最初想法同样有效,并且是手工解决考试练习的一个很好的练习。

登录评论。

标签

世界杯预选赛小组名单社区寻宝

在MATLAB Central中找到宝藏,并发现社区如何帮助您!世界杯预选赛小组名单

开始狩猎!

Baidu
map