【答案】应助回帖

1楼 2011-11-14 22:13:43

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

change0618

铁杆木虫 (著名写手)

方丈大师

计算强帖: 1
应助: 44 (小学生)
金币: 17724.5
红花: 17
帖子: 2413
在线: 546.7小时
虫号: 496517
注册: 2008-01-19
专业: 化学反应工程

【答案】应助回帖

★ ★
dbb627(金币+2): 感谢应助，这个是奇异值问题，你的解法是对的 2011-11-14 23:34:12
dbb627:编辑内容 2011-11-14 23:35

CODE:

function Exercise110

clear all

clc

a = 0;

b = 1;

solinit = bvpinit(linspace(a,b,11),[0 0]);

S = diag([0 -1]);     % 奇异非线性微分方程组y' = S*y/x+f(x,y)的S矩阵,则ODEfun = f(x,y)

options = bvpset('SingularTerm',S);  % 设置奇异性参数S矩阵选项

sol = bvp4c(@ODEfun,@BCfun,solinit,options);

x = 0:0.02:1;

y = deval(sol,x);

plot(x,y(1,:),'ro-')

xlabel('x')

ylabel('y')

function dy=ODEfun(x,y)

dy1=y(2);

dy2=2*y(1)/(1+y(1));

dy=[dy1;dy2];

function bc = BCfun(ya,yb)

bc = [ya(2);yb(1)-1];

[ Last edited by dbb627 on 2011-11-14 at 23:35 ]

2楼2011-11-14 22:45:37

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

dbb627

荣誉版主 (著名写手)

专家经验: +4
计算强帖: 12
应助: 289 (大学生)
贵宾: 0.589
金币: 24640.4
散金: 551
红花: 61
沙发: 1
帖子: 1246
在线: 1794.8小时
虫号: 149791
注册: 2005-12-29
性别: GG
专业: 污染控制化学
管辖: 计算模拟

【答案】应助回帖

★ ★ ★
cenwanglai(金币+3): 谢谢回复！ 2011-11-15 07:11:26
owxiao(金币+30): 2011-11-15 12:27:16

注意，这里换了一个条件，dy/dx（x=0）=0
因为x=0 ${1/x}{dy/dx}$ 无意义
solinit = bvpinit(linspace(0,1,10),[5;1]);
sol = bvp4c(odefun,bcfun,solinit);
xint = linspace(0,1,30);
yint = deval(sol,xint);
hold on; plot(xint,yint(1,

,'rp');
??? Error using ==> bvp4c at 252
Unable to solve the collocation equations -- a
singular Jacobian encountered

CODE:

odefun=@(x,y)[y(2);2*y(1)/(1+y(1))-y(2)/x];

bcfun=@(ya,yb)[ya(2);yb(1)-1];

solinit = bvpinit(linspace(1e-7,1,10),[5;1]);

sol = bvp4c(odefun,bcfun,solinit);

xint = linspace(1e-7,1,30);

yint = deval(sol,xint);

hold on; plot(xint,yint(1,:),'-rp');

结果如图

The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.

3楼2011-11-14 22:56:01

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

刘伟利

木虫 (文坛精英)

应助: 14 (小学生)
金币: 1881.1
红花: 8
帖子: 14472
在线: 767.5小时
虫号: 1339791
注册: 2011-07-07
性别: GG
专业: 化工系统工程

【答案】应助回帖

★ ★
dbb627(金币+1): 有问题，可以名言 2011-11-14 23:38:04
dbb627(金币-3): 这个属于违规复制回答 2011-11-14 23:52:47

function Exercise110
clear all
clc
a = 0;
b = 1;
solinit = bvpinit(linspace(a,b,11),[0 0]);
S = diag([0 -1]);    % 奇异非线性微分方程组y' = S*y/x+f(x,y)的S矩阵,则ODEfun = f(x,y)
options = bvpset('SingularTerm',S);  % 设置奇异性参数S矩阵选项
sol = bvp4c(@ODEfun,@BCfun,solinit,options);
x = 0:0.02:1;
y = deval(sol,x);
plot(x,y(1,,'ro-')
xlabel('x')
ylabel('y')

function dy=ODEfun(x,y)
dy1=y(2);
dy2=2*y(1)/(1+y(1));
dy=[dy1;dy2];

function bc = BCfun(ya,yb)
bc = [ya(2);yb(1)-1];

如果一窍不通,请先看书学习回来再开口
2011-11-14 22:45:37 2
TA的主页  加收听  加好友  只看回帖  发短消息  送红花引用回复编辑通知版主处理
dbb627
实习版主 (小有名气)

计算强帖: 2
听众: 9
红花: 5
贵宾: 0.012
金币: 2942.6
散金: 73
帖子: 231
在线: 194.3小时
虫号: 149791
注册: 2005-12-29
活跃: 2011-11-14
状态: 在线
性别: GG
专业: 污染控制化学
管辖: 计算模拟
【答案】应助回帖

注意，这里换了一个条件，dy/dx（x=0）=0
因为x=0 无意义
solinit = bvpinit(linspace(0,1,10),[5;1]);
sol = bvp4c(odefun,bcfun,solinit);
xint = linspace(0,1,30);
yint = deval(sol,xint);
hold on; plot(xint,yint(1,,'rp');
??? Error using ==> bvp4c at 252
Unable to solve the collocation equations -- a
singular Jacobian encountered
CODE:  [Copy to clipboard]
odefun=@(x,y)[y(2);2*y(1)/(1+y(1))-y(2)/x];
bcfun=@(ya,yb)[ya(2);yb(1)-1];
solinit = bvpinit(linspace(1e-7,1,10),[5;1]);
sol = bvp4c(odefun,bcfun,solinit);
xint = linspace(1e-7,1,30);
yint = deval(sol,xint);
hold on; plot(xint,yint(1,

,'-rp');
结果如图

赞一下(2人)

4楼2011-11-14 23:04:56

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

owxiao

木虫 (小有名气)

应助: 0 (幼儿园)
金币: 2160.5
红花: 1
帖子: 180
在线: 104.3小时
虫号: 1196613
注册: 2011-01-26
性别: GG
专业: 环境工程

引用回帖:

1楼: Originally posted by owxiao at 2011-11-14 22:13:43:
求救。之后还有金币。。。

sol=bvp4c(@ODEfun,@BCfun,solinit,options);
??? Error using ==> feval
Undefined command/function 'BCfun'.

Error in ==> bvp4c at 190
testBC = feval(bc,y(:,1),y(:,end),ExtraArgs{:

这边运行好像一直有错误

5楼2011-11-15 00:18:53

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

owxiao

木虫 (小有名气)

应助: 0 (幼儿园)
金币: 2160.5
红花: 1
帖子: 180
在线: 104.3小时
虫号: 1196613
注册: 2011-01-26
性别: GG
专业: 环境工程

引用回帖:

2楼: Originally posted by change0618 at 2011-11-14 22:45:37:
[code]function Exercise110
clear all
clc
a = 0;
b = 1;
solinit = bvpinit(linspace(a,b,11),[0 0]);
S = diag([0 -1]); % 奇异非线性微分方程组y' = S*y/x+f(x,y)的S矩阵,则ODEfun = f(x,y)
optio ...

6楼2011-11-15 00:19:52

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

dbb627

荣誉版主 (著名写手)

专家经验: +4
计算强帖: 12
应助: 289 (大学生)
贵宾: 0.589
金币: 24640.4
散金: 551
红花: 61
沙发: 1
帖子: 1246
在线: 1794.8小时
虫号: 149791
注册: 2005-12-29
性别: GG
专业: 污染控制化学
管辖: 计算模拟

【答案】应助回帖

★ ★ ★
cenwanglai(金币+3): 谢谢回复！ 2011-11-15 07:11:09

引用回帖:

6楼: Originally posted by owxiao at 2011-11-15 00:19:52:
sol=bvp4c(@ODEfun,@BCfun,solinit,options);
??? Error using ==> feval
Undefined command/function 'BCfun'.

Error in ==> bvp4c at 190
testBC = feval(bc,y(:,1),y(:,end),ExtraArgs{:});

...

是你运行的问题，要将所用代码复制，在editor保存为m文件运行
或者你运行下面代码吧，是等价的

CODE:

odefun=@(x,y)[y(2);2*y(1)/(1+y(1))];

bcfun=@(ya,yb)[ya(2);yb(1)-1];

S = diag([0 -1]);    

options = bvpset('SingularTerm',S);

solinit = bvpinit(linspace(0,1,10),[5;1]);

sol = bvp4c(odefun,bcfun,solinit,options);

x = 0:0.02:1;

y = deval(sol,x);

plot(x,y(1,:),'ro-')

xlabel('x')

ylabel('y')

[ Last edited by dbb627 on 2011-11-15 at 09:18 ]

The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.

7楼2011-11-15 00:40:37

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

dingd

铁杆木虫 (职业作家)

计算强帖: 4
应助: 1641 (讲师)
金币: 15037.3
散金: 101
红花: 234
帖子: 3410
在线: 1223.7小时
虫号: 291104
注册: 2006-10-28

★
dbb627(金币+1): 欢迎交流 2011-12-06 08:48:59

1stOpt求解：

CODE:

Variable x=[1e-10,1],y=[,1], y'=0;

Plot y;

ODEFunction y''=2*y/(1+y)-1/x*y';

8楼2011-12-05 12:06:32

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

mlxqz

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 24
帖子: 6
在线: 17.9小时
虫号: 7600233
注册: 2017-12-10
专业: 机械结构强度学

引用回帖:

2楼: Originally posted by change0618 at 2011-11-14 22:45:37
function Exercise110
clear all
clc
a = 0;
b = 1;
solinit = bvpinit(linspace(a,b,11),);
S = diag(); % 奇异非线性微分方程组y' = S*y/x+f(x,y)的S矩阵,则ODEfun = f(x,y)
options = bvpset('Singul ...

您好，这个S矩阵是如何得出的，谢谢

9楼2018-12-22 15:40:33

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

mlxqz

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 24
帖子: 6
在线: 17.9小时
虫号: 7600233
注册: 2017-12-10
专业: 机械结构强度学

solinit = bvpinit(linspace(0,14,10),[1;0]);
sol=bvp4c(@ODEfun,@BCfun,solinit);
x = linspace(0,14,127);
y = deval(sol,x);
plot(x,y(2,

)

function dydx=ODEfun(x,y)
a=pi*14.701^4/64;
b=pi*(24^4-14.701^4)/64;
c=3.54365802196686;
d=5.94058826800217;
e=2/14.701/pi;
f=14.701/2;
dydx(1,1)=y(2);
dydx(2,1)=(1/a+1/b)/(c+d)*f^2/e*(cos(x/(f*e)))^2*y(1)-(tan(x/(f*e)))/f/e*y(2);

function res=BCfun(ya,yb)
res=[ya(1); yb(1)-1];
哪位大神给看看我这个微分方程的数值解是多少。用bvp4c解出来的总感觉不对