如何用matlab把两个方程整合一个

如何用matlab把两个方程整合一个

这个怎么变成 r与seita的关系式啊 

需要r与θ的关系式 然后有r的数值 求θ

一般这样的微分方程用ode45等函数求解
解得r和t的关系,theta和t的关系
就得到
r和theta的关系了,注意这个关系是两列数值的关系而不是具体的函数形式
对于一些数列中没有的值,可以用插值的办法得到对应的值
由于方程比较复杂,只能用数值解
所以一定要先给出 F m mu v等参数的数值
还有r,theta和r,theta一阶导数的初始值
定一个适合的时间t的取值范围
才能 给出结果
如果你给出这些参数的具体数值,我就可以写出程序
再问： 如果能快点的话 我追加

再答： 初始值呢？ 也就是时间为0时，r和theta以及它们对时间一阶导数的值
没有初始值也解不了的
再问： 这次应该对了 上次写错了

再答： 2.4t是2.4吨的意思，还是2.4乘以时间t？
再问： 吨
再答： mu=4.9x10^23 是什么？怎么这么大？单位是什么？
物理意义是什么？

r的单位是米吗？

由于是数值解，所有数值都不能有问题，所以要问清楚一点
再问： 万有引力常数乘以质量  r单位是千米

再答： 为了得到标准的单位，这里我用国际标准单位制 就是 kg m s制

看上去 7.3477e22kg是月球的质量

mu=GM=6.64e-11 x 7.3477e22 ~ 4.9e12 (m^3/s^2) （这里怎么算出23次方的？）

难怪之前算的结果都很奇怪，等我改好数值再慢慢算
再问： 抱歉啊 这个程序大致几点能出来啊 出来的数值如果拟合成线的话是类似抛物线么
再答： F=7500;mu=4.9e12;m=2.4e3;x0=[1752.013e3    1700      0       0.92465067]';           %r(0)          r'(0)   theta(0)    theta'(0)          fun=@(t,x) [x(2);-mu/x(1)+x(1)*x(4)^2;x(4);-F/m/x(1)-2*x(2)*x(4)/x(1)];ts=[0:0.01:2.8];
[t x]=ode45(fun,ts,x0);r=x(:,1);theta=x(:,3);subplot(2,2,1),plot(t,r);set(gca,'xlim',[ts(1) ts(end)]);ylabel('r(m)');subplot(2,2,3),plot(t,theta);set(gca,'xlim',[ts(1) ts(end)]);ylabel('\theta');xlabel('time(s)');subplot(2,2,[2 4]),plot(theta,r);set(gca,'xlim',[theta(1) theta(end)]);ylabel('r(m)');xlabel('\theta');

以上数据单位都用标准单位，那么时间的单位应该是秒试了一下，在0~2.8秒内，r(theta)的关系类似你所说的抛物线画出来三个曲线表示这0~2.8秒内，r theta t之间的关系
再问： 很清楚 谢谢你了 如果我有什么不懂的问题还可以问你么
再答： 我又看了一下，你的程序初值是有问题的
r0=1752.015e3 m V(0)=1700m/s 是初始半径和初始速度
而这个初始的速度应该是沿着 theta方向的，也就是角方向的
r'的初值应该是0
theta' 的初值是 角速度的初值 等于 V0/r0

我改了一下程序 ，这样结果合理多了，应该是正确的
着陆的时间应该是260秒左右

F=7500;mu=4.9e12;m=2.4e3;
r0=1752.013e3;V0=1700;
x0=[r0 0 0 V0/r0 ]';
%r(0) r'(0) theta(0) theta'(0) 

fun=@(t,x) [x(2);-mu/x(1)^2+x(1)*x(4)^2;x(4);-F/m/x(1)-2*x(2)*x(4)/x(1)];
ts=[0:1:260];

[t x]=ode45(fun,ts,x0);

r=x(:,1);
theta=x(:,3);
rmoon=1738e3;
[xx,yy]=pol2cart(theta,r);
[xx0,yy0]=pol2cart(theta,rmoon);
subplot(2,2,1),plot(t,r);set(gca,'xlim',[ts(1) ts(end)]);
ylabel('r(m)');xlabel('time(s)');
subplot(2,2,3),plot(t,theta);set(gca,'xlim',[ts(1) ts(end)]);
ylabel('\theta');xlabel('time(s)');
subplot(2,2,2),plot(theta,r);set(gca,'xlim',[theta(1) theta(end)]);
ylabel('r(m)');xlabel('\theta');
subplot(2,2,4),
plot(xx0,yy0,'k','linewidth',2);hold on;
plot(xx,yy);hold off;legend('月球表面','轨道');axis equal;
axis([1.6e6 1.6e6+max(yy0) min(yy) max(yy)]);
ylabel('y(m)');xlabel('x(m)');