24小时热门版块排行榜    

查看: 621  |  回复: 4

pizi7880

木虫 (正式写手)

[求助] 数值找根,总感觉对,但总是错。

在求根的时候遇到一个困难:
给定一个2阶倒数微分方程:  z''[t]=z[t]+1/t。在t0时刻初始条件知道,z[t0]=t0,z'[t0]=t0^2。那么我们可以很容易通过NDSOLVE命令得到任意一个时刻 tfin,z[tfin]的值的大小。
而我要做的是,我想找到一个合适t0,使得上面计算得到的z'[tfin]是一个我给定的值 2,即要求,z'[tfin]=2.
下面是我的程序,但是程序总是报错,赐教!
数值找根,总感觉对,但总是错。
无标题.jpg



[ Last edited by pizi7880 on 2013-7-6 at 23:53 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

walk1997

金虫 (著名写手)

vf[t0_?NumericQ,tfin_?NumericQ]
大致是这样的逻辑:数值
具体 自己敲一下
2楼2013-07-07 11:05:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mshwangg

至尊木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
pizi7880: 金币+8 2013-07-07 16:40:23
xzhdty: 金币+1, 感谢参与 2013-07-07 21:52:40
NDSolve求解出来的是z[t],直接替换z'[t]是不可以的。需要替换出来z[t],自己再做微分。
其实这是一个打靶的问题
tfirst = 5;
tfin = 10;
eq[t0_] := {z''[t] == z[t] + 1/t, z[t0] == t0, z'[t0] == t0^2};
f[t0_] := NDSolve[eq[t0], z[t], {t, t0, tfin}]
y[t0_?NumericQ] := First[D[z[t] /. f[t0], t] /. t -> tfin]
r = FindRoot[y[t0] - 2 == 0, {t0, tfirst, tfin}]
Plot[Evaluate[z[t] /. f[t0 /. r]], {t, t0 /. r, tfin}]
Plot[Evaluate[D[z[t] /. f[t0 /. r], t]], {t, t0 /. r, tfin}]

计算出来{t0 -> 9.90643},有一定误差
3楼2013-07-07 12:55:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

walk1997

金虫 (著名写手)

★ ★ ★ ★ ★ ★
pizi7880(xzhdty代发): 金币+6 2013-07-07 21:52:58
内容已删除
4楼2013-07-07 13:14:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mshwangg

至尊木虫 (正式写手)

【答案】应助回帖

引用回帖:
4楼: Originally posted by walk1997 at 2013-07-07 13:14:55
NDSolve出来的是可以直接微分的 把自变量换下
关键是NumericQ的使用
对的 是个打靶的问题
不过你上面这个结果 好像不太对  其实已经发散了(FindRoot 没有收敛)
LZ的这个问题比较特殊 方程其实有解析解.....
...

你说的有道理。
FindRoot 100个循环已经不够用了。所以存在误差的。
show一下解析解吧,验证一下
{{z[t] ->
   1/2 E^(-t -
     t0) (E^(2 t) t0 + E^(2 t0) t0 + E^(2 t) t0^2 - E^(2 t0) t0^2 +
      E^(2 t + t0) ExpIntegralEi[-t] - E^t0 ExpIntegralEi[t] -
      E^(2 t + t0) ExpIntegralEi[-t0] + E^t0 ExpIntegralEi[t0])}}
tfin=1.46时,z'[t0]关于t0作图
数值找根,总感觉对,但总是错。-1
可看到有两个t0是满足要求的,walk的结果是合理的。
而tfin=10的时候,z'[t0]与2那条直线无交点,所以z'[t0]=2无解
5楼2013-07-07 16:49:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 pizi7880 的主题更新
信息提示
请填处理意见