24小时热门版块排行榜    

查看: 939  |  回复: 7
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

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编程绘图

» 猜你喜欢

爆炸力学
已阅   回复此楼   关注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的回帖
查看全部 8 个回答

月只蓝

主管区长 (职业作家)

【答案】应助回帖

感谢参与,应助指数 +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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 280求调剂 +7 咕噜晓晓 2026-03-18 8/400 2026-03-21 01:27 by JourneyLucky
[考研] 295求调剂 +4 一志愿京区211 2026-03-18 6/300 2026-03-20 23:41 by JourneyLucky
[考研] 304求调剂 +7 司空. 2026-03-18 7/350 2026-03-20 23:08 by JourneyLucky
[考研] 一志愿武汉理工材料工程专硕调剂 +9 Doleres 2026-03-19 9/450 2026-03-20 22:36 by JourneyLucky
[考研] 290求调剂 +7 ^O^乜 2026-03-19 7/350 2026-03-20 21:43 by JourneyLucky
[考研] 北科281学硕材料求调剂 +5 tcxiaoxx 2026-03-20 5/250 2026-03-20 21:35 by laoshidan
[考研] 一志愿西南交通 专硕 材料355 本科双非 求调剂 +5 西南交通专材355 2026-03-19 5/250 2026-03-20 21:10 by JourneyLucky
[考研] 广西大学家禽遗传育种课题组2026年硕士招生(接收计算机专业调剂) +3 123阿标 2026-03-17 3/150 2026-03-20 15:58 by 飞行琦
[考研] 298-一志愿中国农业大学-求调剂 +9 手机用户 2026-03-17 9/450 2026-03-20 14:24 by 无懈可击111
[考研] 能源材料化学课题组招收硕士研究生8-10名 +5 脱颖而出 2026-03-16 14/700 2026-03-20 09:30 by kkcoco25
[考研] 一志愿中国海洋大学,生物学,301分,求调剂 +5 1孙悟空 2026-03-17 6/300 2026-03-19 23:46 by zcl123
[考研] 081700化工学硕调剂 +3 【1】 2026-03-16 3/150 2026-03-19 23:40 by edmund7
[考研] 0703化学调剂 +10 妮妮ninicgb 2026-03-15 14/700 2026-03-19 22:59 by 学员8dgXkO
[考研] 085601材料工程专硕求调剂 +10 慕寒mio 2026-03-16 10/500 2026-03-19 15:26 by 丁丁*
[考研] 收复试调剂生 +4 雨后秋荷 2026-03-18 4/200 2026-03-18 14:16 by elevennnne
[考研] 301求调剂 +4 A_JiXing 2026-03-16 4/200 2026-03-17 17:32 by ruiyingmiao
[考研] 308求调剂 +4 是Lupa啊 2026-03-16 4/200 2026-03-17 17:12 by ruiyingmiao
[考研] 有没有道铁/土木的想调剂南林,给自己招师弟中~ +3 TqlXswl 2026-03-16 7/350 2026-03-17 15:23 by TqlXswl
[考研] 0854控制工程 359求调剂 可跨专业 +3 626776879 2026-03-14 9/450 2026-03-16 17:42 by 626776879
[考研] 本科南京大学一志愿川大药学327 +3 麦田耕者 2026-03-14 3/150 2026-03-14 20:04 by 外星文明
信息提示
请填处理意见