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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » MATLAB/Mathematica » 这个程序是不是哪里有问题？初始值跟实验值是一样的，但是模拟结果不对
81/1返回列表
查看: 929  |  回复: 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返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 学校已经提交到NSFC，还能修改吗？ 40+4 babangida 2026-03-19 8/400 2026-03-20 15:58 by babero
[考研] 296求调剂 +4 www_q 2026-03-18 8/400 2026-03-20 12:12 by 学员8dgXkO
[考研] 一志愿西南交通 专硕 材料355 本科双非 求调剂 +4 西南交通专材355 2026-03-19 4/200 2026-03-20 11:39 by 花开富贵幸福人�
[考研] 能源材料化学课题组招收硕士研究生8-10名 +5 脱颖而出 2026-03-16 14/700 2026-03-20 09:30 by kkcoco25
[考博] 申博26年 +3 八6八68 2026-03-19 3/150 2026-03-19 19:43 by nxgogo
[考研] 0703化学调剂 ，六级已过，有科研经历 +12 曦熙兮 2026-03-15 12/600 2026-03-19 19:42 by maocaozhuxi
[考博] 东华理工大学化材专业26届硕士博士申请 +8 zlingli 2026-03-13 8/400 2026-03-19 16:32 by 轻松不少随
[考研] 材料与化工求调剂 +7 为学666 2026-03-16 7/350 2026-03-19 14:48 by 尽舜尧1
[考研] 085600材料与化工调剂 324分 +10 llllkkkhh 2026-03-18 12/600 2026-03-19 14:33 by llllkkkhh
[考研] 化学求调剂 +3 临泽境llllll 2026-03-17 4/200 2026-03-19 13:59 by houyaoxu
[考研] 一志愿中海洋材料工程专硕330分求调剂 +7 小材化本科 2026-03-18 7/350 2026-03-19 10:46 by Linda Hu
[考研] 材料工程专硕调剂 +5 204818@lcx 2026-03-17 6/300 2026-03-18 22:55 by 204818@lcx
[考研] 344求调剂 +6 knight344 2026-03-16 7/350 2026-03-18 20:13 by walc
[考研] 302求调剂 +10 呼呼呼。。。。 2026-03-17 10/500 2026-03-18 12:45 by Linda Hu
[考研] 0703化学调剂 +3 妮妮ninicgb 2026-03-17 3/150 2026-03-18 10:29 by macy2011
[考研] 环境工程调剂 +8 大可digkids 2026-03-16 8/400 2026-03-18 09:36 by zhukairuo
[基金申请] 今年的国基金是打分制吗？ 50+3 zhanghaozhu 2026-03-14 3/150 2026-03-16 17:07 by 北京莱茵润色
[考研] 26考研一志愿中国石油大学(华东)305分求调剂 +3 嘉年新程 2026-03-15 3/150 2026-03-15 13:58 by 哈哈哈哈嘿嘿嘿
[考研] 297一志愿上交085600求调剂 +5 指尖八千里 2026-03-14 5/250 2026-03-14 17:26 by a不易
[考研] 297求调剂 +4 学海漂泊 2026-03-13 4/200 2026-03-14 11:51 by 热情沙漠
信息提示
关闭
请填处理意见
关闭