吴老师，您好。

首先，谢谢吴老师为大家提供一本这么好的教材，你的书太给力了。

我最近在求解一个三阶常微分方程边值问题，您的书帮助很大。您书的第十章有关于边值问题的例子，跟我求解的问题很相似，所以我是直接在例子10.6-3的程序代码的基础上做的。结果出来了，但是结果跟所求问题在物理上不符合，自己研究了很久没看找出问题出在哪，麻烦吴老师指导一下。先谢谢您了。

我要求的常微分方程是：

方程1.JPG (21.58 KB, 下载次数: 29)

常微分方程

2011-11-17 21:53 上传

F是eta的函数，F‘,F’‘分别是eta的一阶，二阶导数。使用bvp4c求解，因为有四个边界，所以我参照例子10.6.3引入一个参数lambda

把方程中常数项设成lambda=1.8348

结果(图1)存在的问题是：

运行结果F’的变化趋势与想要求的结果不一致。结果图中F’的曲线是向下凹的。我想要的结果是向上凸的图。

因为要用F’计算另外的参数u,y，

从图1运行结果得到的u与y的关系图(图2)与想要的结果(图3)不一致。

方程2_u,v.JPG (4.49 KB, 下载次数: 9)

u，y与eta，F‘的关系

2011-11-17 21:51 上传

书求助1.jpg (22.93 KB, 下载次数: 10)

常微分方程求解结果

2011-11-17 21:51 上传

图1 常微分方程求解结果

书求助2.jpg (18.07 KB, 下载次数: 10)

由常微分方程得到的u与y的关系图

2011-11-17 21:51 上传

图2 由常微分方程得到的u与y的关系图

书求助3.JPG (17.49 KB, 下载次数: 5)

想要得到的u与y的关系图

2011-11-17 21:51 上传

图3 想要得到的u与y的关系图

我用的程序代码：

function example10_6_3

n=47.4256;

X = linspace(0,n,1000);%为bvpinit生成初始化网络准备

lambda = 1.8348;%对lambda的猜测值

function yinit = mat4init(x)

%对y初值的估计函数,由于y2(0) = 1,y2(pi) = 0;所以挑选一个满足上述条件的函数

%这里选择的是-k*sin(4*x)来作为对y2(t)的估计，其原函数y1(t)满足y1(0)=1，因此

%k = 4，所以y1(t)的估计为cos(4*x),y2(t)的估计为-4*sin(4*x)

yinit = [-(8.8920998e-4)*x^3+(0.063256975)*x^2;-(8.8920998e-4)*3*x^2+(0.063256975)*2*x;-(8.8920998e-4)*6*x+(0.063256975)*2];

end

function res = mat4bc(ya,yb,lambda)

%lambda作为未知参数不能省略

%边界条件;边界条件为y1(0) = 1,y2(0) = 2,y2(pi) = 0; 0,0,pi分别对应ya,ya,yb

%而y1,y2,y2分别是第一第二第二个状态变量。因此得到如下res表达式

res = [ya(1)

ya(2)

yb(1)-47.4256

yb(2)];

end

function dydx = mat4ode(x,y,lambda)

%微分方程函数

dydx = [ y(2);y(3);-y(1)*y(3)/2+y(2)^2/2+0.8348*y(1)+lambda];

end

solinit = bvpinit(X,@mat4init,lambda);%由bvpinit生成的初始化网格

sol = bvp4c(@mat4ode,@mat4bc,solinit);%调用bvp4c求解,也可以换成bvp5c

tint = linspace(0,n);

Stint = deval(sol,tint);%根据得到的sol利用deval函数求出[0,4]区间内更多其他的解

%画图展示

figure;

plot(tint,Stint(1,:),'k-','linewidth',2);

hold on

plot(tint,Stint(2,:),'k:','linewidth',2);

hold on

plot(tint,Stint(3,:),'k--','linewidth',2); %这里是(3, : ) 不知道怎么回事网页上显示笑脸

L = legend('{\itF}(\eta)','{\itdF}(\eta)','{\itddF}(\eta)','Location','best');

set(L,'fontname','Times New Roman');

xlabel('\itt','fontsize',16);ylabel('方程的解');

-Stint(2,1)*Stint(3,1)/2+Stint(2,1)^2/2+0.8347708*Stint(1,1)+lambda

(Stint(3,2)-Stint(3,1))/(n/5000)

end

[本帖最后由 吴鹏 于 2011-11-21 14:08 编辑]