应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » 这个程序是不是哪里有问题？初始值跟实验值是一样的，但是模拟结果不对
81/1返回列表
查看: 1188  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

zhaoqian59

捐助贵宾 (小有名气)

[求助] 这个程序是不是哪里有问题？初始值跟实验值是一样的，但是模拟结果不对 已有1人参与

f(x)
pl=1000;
phyd=110000;
P0=29.39E9;
a0=0.0341;
v=1.5;

x1=x(1);
x2=x(2);
pg=P0*(a0/x1)^(3*v);
y=[x2 (1/pl*(pg-phyd)-3/2*x2^2)/x1]';

Ts = 0.001;    % 仿真步长
Tn = 1;       % 仿真终止时间
t  = 0:Ts:Tn;   % 仿真时间范围
N  = length(t);
x=[0.034 0]'; % 系统初值：x=[a a的导数]
%% Variable Declaration and Initialization
a=0.034*ones(1,N);
a_dt=zeros(1,N);
%x = [a a_d]'; % intial states
%% Simulation Iteration
for n = 1:N
    a(n)=x(1);
    a_dt(n)=x(2);
    k1 = fx(x);
    k2 = fx(x+Ts/2*k1);
    k3 = fx(x+Ts/2*k2);
    k4 = fx( x+Ts *k3);
    x = x + Ts/6 * ( k1 + 2*k2 + 2*k3 + k4 );
end
%% Plot Response
plot(t,a)
xlabel('t/(s)')
ylabel('a')
是一个二阶非线性常微分方程
回复此楼

» 收录本帖的淘帖专辑推荐

matlab编程绘图

» 猜你喜欢
爆炸力学
1楼 2015-08-11 09:29:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

感谢参与，应助指数 +1
楼主自己编写的四阶龙格库塔来求解常微分方程？
能给出方程具体形式么，初值，以及各常数？
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
2楼2015-08-11 10:37:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoqian59

捐助贵宾 (小有名气)
引用回帖:
2楼: Originally posted by 月只蓝 at 2015-08-11 10:37:11
楼主自己编写的四阶龙格库塔来求解常微分方程？
能给出方程具体形式么，初值，以及各常数？

aa''+1.5a'=1/1000*(P-Pwater)
P=P0*(a0/a)^v
a是关于t的函数，Pwater=1.1*100000，P0=29E9，a0=0.0341，v=3.谢谢你啦~
赞一下 回复此楼
爆炸力学
3楼2015-08-11 15:03:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)
引用回帖:
3楼: Originally posted by zhaoqian59 at 2015-08-11 15:03:26
aa''+1.5a'=1/1000*(P-Pwater)
P=P0*(a0/a)^v
a是关于t的函数，Pwater=1.1*100000，P0=29E9，a0=0.0341，v=3.谢谢你啦~...

t=0时，对应a a'的初值呢？
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
4楼2015-08-11 15:08:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoqian59

捐助贵宾 (小有名气)
引用回帖:
4楼: Originally posted by 月只蓝 at 2015-08-11 15:08:12
t=0时，对应a a'的初值呢？...

t=0时，a(t)=0.0341,a'=0
赞一下 回复此楼
爆炸力学
5楼2015-08-11 15:14:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

引用回帖:
3楼: Originally posted by zhaoqian59 at 2015-08-11 15:03:26
aa''+1.5a'=1/1000*(P-Pwater)
P=P0*(a0/a)^v
a是关于t的函数，Pwater=1.1*100000，P0=29E9，a0=0.0341，v=3.谢谢你啦~...

记 a=u1，a'=u2；
原二阶方程可将阶为方程组：
u1'=u2;
u2'=1/u1*  (  (P-Pwater)/1000 - 1.5*u2     );
其中，P=P0*(a0/u1)^v;

MATLAB代码如下：
CODE:
function solve_odes_two_oder
clear all;clc
u0=[0.0341 0];
tspan=linspace(0,0.1,100);

[t u]=ode45(@odefun,tspan,u0);

figure(1)
plot(t,u(:,1),'bo--',t,u(:,2),'r-*'),legend('a',' da/dt ')




function f=odefun(t,u)
Pwater=1.1*100000;
P0=29E9;
a0=0.0341;
v=3;
P=P0*(a0/u(1))^v;

f(1)=u(2);
f(2)=1/u(1)*   (  (P-Pwater)/1000 - 1.5*u(2)     );
f=f';

结果中，a和a'迅速增大，即便在0.1的时间之内，见附图1。
这个程序是不是哪里有问题？初始值跟实验值是一样的，但是模拟结果不对
附图1.png
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
6楼2015-08-11 15:47:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhaoqian59

捐助贵宾 (小有名气)
引用回帖:
6楼: Originally posted by 月只蓝 at 2015-08-11 15:47:29
记 a=u1，a'=u2；
原二阶方程可将阶为方程组：
u1'=u2;
u2'=1/u1*  (  (P-Pwater)/1000 - 1.5*u2     );
其中，P=P0*(a0/u1)^v;

MATLAB代码如下：

function solve_odes_two_oder
clear all;clc
u0=;
...

复制过去之后运行不出来.....
赞一下 回复此楼
爆炸力学
7楼2015-08-12 09:40:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

月只蓝

主管区长 (职业作家)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
zhaoqian59: 金币+50 2015-08-12 12:40:35
引用回帖:
7楼: Originally posted by zhaoqian59 at 2015-08-12 09:40:23
复制过去之后运行不出来........

代码完整复制到一个新建的m文件，运行即可，不要在主程序窗口运行。

[ 发自小木虫客户端 ]
赞一下 回复此楼
MATLAB、MS小问题、普通问题请发帖求助！时间精力有限，恕不接受无偿私信求助。
8楼2015-08-12 10:17:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zhaoqian59 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[找工作] 售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急 +4 ZYXdhzDAy9ZX 2026-09-28 4/200 2026-09-29 16:54 by yCO1Ll7aHtsw
[考研] 售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急 +3 ZYXdhzDAy9ZX 2026-09-28 4/200 2026-09-29 16:37 by etmYJ6d2rquH
[博后之家] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 ZYXdhzDAy9ZX 2026-09-28 3/150 2026-09-29 16:30 by etmYJ6d2rquH
[硕博家园] 售SCI文章，我:8O.5.5.1O.54,科目全,可十急 +3 ZYXdhzDAy9ZX 2026-09-28 3/150 2026-09-29 16:03 by etmYJ6d2rquH
[考博] 售一区SCI文章T0P，我:8O.551.O54,科目全,可十急 +3 tqUTOxClQUMF 2026-09-28 3/150 2026-09-29 15:46 by etmYJ6d2rquH
[论文投稿] 售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急 +4 tqUTOxClQUMF 2026-09-28 4/200 2026-09-29 15:15 by etmYJ6d2rquH
[硕博家园] 售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急 +3 GDBe8tDZqE8z 2026-09-28 3/150 2026-09-29 14:30 by etmYJ6d2rquH
[考研] 售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急 +3 CfXuS1rDhLYN 2026-09-28 3/150 2026-09-29 14:24 by etmYJ6d2rquH
[找工作] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +6 IDs3scOF0tjC 2026-09-28 6/300 2026-09-29 14:13 by etmYJ6d2rquH
[公派出国] 售一区SCI文章T0P，我:8O.551.O54,科目全,可十急 +4 IDs3scOF0tjC 2026-09-28 5/250 2026-09-29 13:50 by etmYJ6d2rquH
[教师之家] 某top大学教授说“能够在市场中兑现的能力才是真能力”无比同意！ +8 zju2000 2026-09-26 9/450 2026-09-29 13:47 by ahsoarli
[教师之家] 售一区SCI文章T0P，我:8O.551.O54,科目全,可十急 +4 GDBe8tDZqE8z 2026-09-28 4/200 2026-09-29 09:49 by JzYBbHIXWSrW
[找工作] 售SCI一区文章，我:8.O.551.O.5.4,科目全,可伽急 +3 ZYXdhzDAy9ZX 2026-09-28 3/150 2026-09-29 06:28 by ZvyPK8n6Nfki
[硕博家园] 售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急 +3 tqUTOxClQUMF 2026-09-28 3/150 2026-09-28 23:46 by b8fCuTHqckEt
[考研] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 GDBe8tDZqE8z 2026-09-28 3/150 2026-09-28 23:12 by b8fCuTHqckEt
[公派出国] 售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急 +4 IDs3scOF0tjC 2026-09-28 4/200 2026-09-28 22:55 by ez6fGg9abYaj
[硕博家园] 售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急 +3 CfXuS1rDhLYN 2026-09-28 4/200 2026-09-28 22:55 by ez6fGg9abYaj
[有机交流] 同一个分子，一条来自文献，一条来自AI——不告诉你答案，你会选哪条？ +3 tianxiaoxian 2026-09-28 6/300 2026-09-28 22:44 by tianxiaoxian
[考博] 售SCI文章，我:8O.5.5.1O.54,科目全,可十急 +3 IDs3scOF0tjC 2026-09-28 6/300 2026-09-28 22:15 by ez6fGg9abYaj
[论文投稿] 售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急 +3 CfXuS1rDhLYN 2026-09-28 3/150 2026-09-28 19:08 by 9lS3ad5oOymn
信息提示
关闭
请填处理意见
关闭