版块导航: 正在加载中...

客户端APP下载

登录注册

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

24小时热门版块排行榜

gzqdyouxia

金虫 (著名写手)

应助: 1 (幼儿园)
贵宾: 0.393
金币: 928.9
帖子: 1600
在线: 310.6小时
虫号: 802190

[交流] 【求助】matlab求解非线性方程组，并画图处理。要求y,z是实数解！

y=z+(7999022433351517/2305843009213693952*z+5509600814111415/72057594037927936*z^3)/(1-a)+1/2*b；
z=y+(-4534249478333089/576460752303423488*y-3611543935826555/576460752303423488*y^3+17523/50000*y^5)/a-1/2*b
注：y，z是因变量，a，b是自变量，其中a的变化范围是（0,1），b的变化范围是（0,0.1）。要求分别画出以a，b为自变量轴，y，z为因变量轴的三维图（即y-a-b，z-a-b）。小虫还有一个问题就是让b为固定值，如b=0.01，怎么画出y-a，z-a（其中a的变化范围是（0,1））的二维图？

[ Last edited by gzqdyouxia on 2010-12-27 at 19:37 ]

» 猜你喜欢

售SCI一区T0P文章，我:8O.55.1.O.5.4,科目齐全,可+急已经有6人回复
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急已经有4人回复
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急已经有7人回复
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急已经有6人回复
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急已经有6人回复
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急已经有6人回复
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急已经有5人回复
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急已经有4人回复
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急已经有6人回复
售SCI文章，我:8O.5.5.1O.54,科目全,可十急已经有4人回复

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

MATLAB用于三元回归参数已经有3人回复
有没有什么办法能控制uds求解的标量在整个求解域上的范围啊已经有6人回复
高斯09w做过渡金属模拟，请教出现的L103错误如何解决，断电中断后这样处理对不对已经有3人回复
matlab中二维图作图已经有8人回复
关于污水处理的问题已经有32人回复
如何使用matlab求解非线性方程组的所有整数解？已经有9人回复
Matlab求解偏微分方程组已经有22人回复
matlab求解非线性方程组已经有16人回复
matlab数值求解边界条件微分方程组已经有7人回复
求matlab求解一个方程组的问题已经有6人回复
求高人指点用matlab求解非线性方程组，解决了追加100金币；已经有11人回复
matlab的fsove 命令求解非线性方程组已经有6人回复
【求助】向各位大侠求助matlab求解微分方程组遇到的一个问题已经有21人回复
【求助】matlab怎么求解偏微分方程组啊，先谢谢了已经有13人回复
【求助】请教虫友关于matlab差分方程的求解和画图程序已经有4人回复
【求助】求教matlab解非线性方程组已经有9人回复

» 抢金币啦！回帖就可以得到:

查看全部散金贴

1楼 2010-12-27 16:51:57

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

xiegangmai

版主 (职业作家)

★ ★ ★
gzqdyouxia(金币+50):期待您进一步指教 2010-12-28 00:21:14
robert2020(金币+3):辛苦了！ 2010-12-29 09:34:34

a 和 b 之间有什么关系吗？

如果 a 和 b 一一对应，可以用fsolve求解，大概如下：

函数：

CODE:

function F = gzqdyouxia( x, a, b )

F = [ - x( 1 )  + x( 2 ) + ( 7999022433351517 / 2305843009213693952 * x( 2 ) + 5509600814111415 / 72057594037927936 * x( 2 ) ^ 3 ) / ( 1 - a ) + 1 / 2 * b;

    - x( 2 ) + x( 1 ) + ( - 4534249478333089 / 576460752303423488 * x( 1 ) - 3611543935826555 / 576460752303423488 * x( 1 ) ^ 3 + 17523 / 50000 * x( 1 ) ^ 5 ) / a - 1 / 2 * b ];

求解：

CODE:

clear

clc

N = 1000;

a = 0 : 1 / N : 1;

b = 0 : 0.1 / N : 0.1;

y = zeros( 1, N );

z = zeros( 1, N );

F0 = [ 0; 0 ];

options=optimset('Display','off');

for i = 2 : N    % a 不取0， 也可以从0开始

    F0 = fsolve( @( x )gzqdyouxia( x, a( i ), b( i ) ), F0, options );

    y( i ) = F0( 1 );

    z( i ) = F0( 2 );

end

plot3( a( 2: end ), b( 2: end ), y, 'r' )

hold on

plot3( a( 2: end ), b( 2: end ), z, 'b' )

xlabel( 'a' )

ylabel( 'b' )

legend( 'y', 'z' )

要画二维曲线，b 为固定值，a在区间递增，求解方程，画图。

如果a 、b没有对应关系，求解出的应该是三维的曲面，大概如下：

CODE:

clear

clc

N = 100;

[ a, b ] = meshgrid( 0 : 1 / N : 1, 0 : 0.1 / N : 0.1);

[ m, n ] = size( a );

y = zeros( m, n );

z = zeros( m, n );

F0 = [ 0; 0 ];

options=optimset('Display','off');

for i = 1 : m

    for j = 1 : n

        F0 = fsolve( @( x )gzqdyouxia( x, a( i, j ), b( i, j ) ), F0, options );

        y( i, j ) = F0( 1 );

        z( i, j ) = F0( 2 );

    end

end

surf( a, b, y )

hold on

surf( a, b, z )

xlabel( 'a' )

ylabel( 'b' )

在options选项中，可以设置求解算法，具体可参考帮助中的fsolve，我试了试用options=optimset('Display','off', 'ScaleProblem', 'Jacobian' );结果略有不同。

[ Last edited by xiegangmai on 2010-12-27 at 20:20 ]

赞一下(1人)

2楼2010-12-27 19:40:23

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

gzqdyouxia

金虫 (著名写手)

应助: 1 (幼儿园)
贵宾: 0.393
金币: 928.9
帖子: 1600
在线: 310.6小时
虫号: 802190

引用回帖:

Originally posted by xiegangmai at 2010-12-27 19:40:23:
a 和 b 之间有什么关系吗？

如果 a 和 b 一一对应，可以用fsolve求解，大概如下：

函数：
[code]function F = gzqdyouxia( x, a, b )

F = [ - x( 1 ) + x( 2 ) + ( 7999022433351517 / 23058430092136 ...

谢谢你的回答，不过问题我没有解决，a，b都是自变量，y，z是关于a和b的二元函数，我要求的是实数对应关系，您能不能棒我再看一下？

3楼2010-12-27 22:18:05

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

xiegangmai

版主 (职业作家)

★ ★
wuming524(金币+2):谢谢交流，辛苦了 2010-12-28 01:23:32

引用回帖:

Originally posted by gzqdyouxia at 2010-12-27 22:18:05:

谢谢你的回答，不过问题我没有解决，a，b都是自变量，y，z是关于a和b的二元函数，我要求的是实数对应关系，您能不能棒我再看一下？

plot(a,y)
plot(b,z)
就可以了啊

赞一下(1人)

4楼2010-12-27 22:56:06

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

gzqdyouxia

金虫 (著名写手)

应助: 1 (幼儿园)
贵宾: 0.393
金币: 928.9
帖子: 1600
在线: 310.6小时
虫号: 802190

引用回帖:

Originally posted by xiegangmai at 2010-12-27 22:56:06:

plot(a,y)
plot(b,z)
就可以了啊

二维的怎么画？我是个初学者，您能不能教我一下，b=0.01

5楼2010-12-27 23:08:04

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

小木虫:)

荣誉版主 (著名写手)

应助: 0 (幼儿园)
贵宾: 2.542
金币: 9239.2
帖子: 1211
在线: 343.9小时
虫号: 800712

★ ★
robert2020(金币+2):多谢应助！辛苦了！ 2010-12-29 09:34:53

引用回帖:

Originally posted by gzqdyouxia at 2010-12-27 23:08:04:

二维的怎么画？我是个初学者，您能不能教我一下，b=0.01

要画a-y图就要固定b和z这里我们假设z=1.2

程序为

clear
clc
b=0.01;
z=1.2;
for i=1:1:100
a(i)=(i-1)*0.01;
y(i)=z+(7999022433351517/2305843009213693952*z+5509600814111415/72057594037927936*z^3)/(1-a(i))+1/2*b;
end
plot(a,y)

注意这里不能使得a=1，否则分母为零画不出图来

要画a-z图就要固定b和y这里我们假设y=1.5

程序为

clear
clc
b=0.01;
y=1.5;
for i=1:1:100
a(i)=i*0.01;
z(i)=y+(-4534249478333089/576460752303423488*y-3611543935826555/576460752303423488*y^3+17523/50000*y^5)/a(i)-1/2*b;
end
plot(a,z)

注意这里不能使得a=0，否则同样无法画图

赞一下(1人)

6楼2010-12-29 03:08:33

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

云淡天高

木虫 (小有名气)

应助: 0 (幼儿园)
金币: 3146.7
帖子: 249
在线: 23.5小时
虫号: 452333

你可以学一下matlab

7楼2010-12-29 08:12:36

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

qwalter

捐助贵宾 (小有名气)

应助: 0 (幼儿园)
金币: 806.9
帖子: 298
在线: 207.2小时
虫号: 1152501

学习下。

8楼2010-12-29 08:21:08

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

blue.wh

铁杆木虫 (知名作家)

应助: 972 (博后)
金币: 8637.6
帖子: 7822
在线: 3306.3小时
虫号: 551495

引用回帖:

Originally posted by xiegangmai at 2010-12-27 19:40:23:
a 和 b 之间有什么关系吗？

如果 a 和 b 一一对应，可以用fsolve求解，大概如下：

函数：
[code]function F = gzqdyouxia( x, a, b )

F = [ - x( 1 ) + x( 2 ) + ( 7999022433351517 / 23058430092136 ...

为什么a不从0开始算起，比如：a = 1.0000e-010 : 1 / N : 1;会出现这个结果：
??? Error using ==> plot
Vectors must be the same lengths.

Error in ==> Untitled at 19
plot( a( 2: end ), y, 'r' )
谢谢~

[ Last edited by blue.wh on 2010-12-29 at 11:24 ]

9楼2010-12-29 11:20:35

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

abingchem

木虫 (著名写手)

应助: 2 (幼儿园)
金币: 4871
帖子: 1672
在线: 117.6小时
虫号: 206219

引用回帖:

Originally posted by xiegangmai at 2010-12-27 19:40:23:
a 和 b 之间有什么关系吗？

如果 a 和 b 一一对应，可以用fsolve求解，大概如下：

函数：
[code]function F = gzqdyouxia( x, a, b )

F = [ - x( 1 ) + x( 2 ) + ( 7999022433351517 / 23058430092136 ...

为什么结果都是零？

10楼2010-12-30 14:12:58

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

blue.wh

铁杆木虫 (知名作家)

应助: 972 (博后)
金币: 8637.6
帖子: 7822
在线: 3306.3小时
虫号: 551495

引用回帖:

Originally posted by xiegangmai at 2010-12-27 19:40:23:
a 和 b 之间有什么关系吗？

如果 a 和 b 一一对应，可以用fsolve求解，大概如下：

函数：
[code]function F = gzqdyouxia( x, a, b )

F = [ - x( 1 ) + x( 2 ) + ( 7999022433351517 / 23058430092136 ...

你用的方法有问题啊，得到结果不是正解！！！！

11楼2011-01-06 20:51:41

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

xiegangmai

版主 (职业作家)

引用回帖:

Originally posted by abingchem at 2010-12-30 14:12:58:

为什么结果都是零？

这和选取的初值有关系。

12楼2011-01-09 16:20:00

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

xiegangmai

版主 (职业作家)

★
wuming524(金币+1):辛苦了 2011-01-09 17:14:27

引用回帖:

Originally posted by blue.wh at 2011-01-06 20:51:41:

你用的方法有问题啊，得到结果不是正解！！！！

不好意思，未能善始善终，呵呵。

不知道这是不是正确的：

CODE:

clear

clc

N = 100;

[ a, b ] = meshgrid( 0 : 1 / N : 1, 0 : 0.1 / N : 0.1 );

[ m, n ] = size( a );

y = zeros( m, n );

z = zeros( m, n );

F0 = [ 0.1; 0.1 ];

options=optimset('Display','off');

for i = 1 : m 

    for j = 1 : n

        F0 = fsolve( @( x )gzqdyouxia( x, a( i, j ), b( i, j ) ), F0, options );

        y( i, j ) = F0( 1 );

        z( i, j ) = F0( 2 );

    end

end

figure( 'name', '三维图')

subplot( 1, 2, 1 )

surf( a, b, y )

title( 'y-a-b' )

subplot( 1, 2, 2 )

surf( a, b, z )

title( 'z-a-b' )

运行过程中有些警告，和方程及初值选取有关。得到的图为：

赞一下(1人)

13楼2011-01-09 16:26:17

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

blue.wh

铁杆木虫 (知名作家)

应助: 972 (博后)
金币: 8637.6
帖子: 7822
在线: 3306.3小时
虫号: 551495

引用回帖:

Originally posted by xiegangmai at 2011-01-09 16:26:17:

不好意思，未能善始善终，呵呵。

不知道这是不是正确的：

[code]clear
clc

N = 100;
[ a, b ] = meshgrid( 0 : 1 / N : 1, 0 : 0.1 / N : 0.1 );
[ m, n ] = size( a );

y = zeros( m, n );
z ...

不对啊，初值改了也不对。。。。。

14楼2011-01-10 10:06:01

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

xiegangmai

版主 (职业作家)

引用回帖:

Originally posted by blue.wh at 2011-01-10 10:06:01:

不对啊，初值改了也不对。。。。。

你的方程组是高阶的非线性方程组，刚度大，用 fsolve求解是有问题的。
方程组应该有多组解，我再用别的解法试试。

关注中。

15楼2011-01-10 10:40:27

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

blue.wh

铁杆木虫 (知名作家)

应助: 972 (博后)
金币: 8637.6
帖子: 7822
在线: 3306.3小时
虫号: 551495

引用回帖:

Originally posted by xiegangmai at 2011-01-10 10:40:27:

你的方程组是高阶的非线性方程组，刚度大，用 fsolve求解是有问题的。
方程组应该有多组解，我再用别的解法试试。

关注中。

是啊，我算了一下二维的比如b=0.01，在a等于0附近，y，z分别对应三组初值，即使用这个三组初值，运行程序，画出图形，只有一组是正解，其余两组都不是。

16楼2011-01-10 10:43:03

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

xiegangmai

版主 (职业作家)

引用回帖:

Originally posted by blue.wh at 2011-01-10 10:43:03:

是啊，我算了一下二维的比如b=0.01，在a等于0附近，y，z分别对应三组初值，即使用这个三组初值，运行程序，画出图形，只有一组是正解，其余两组都不是。

改了一下，计算了一下误差，源码如下：

CODE:

clear

clc

eps = 0.0001;

N = 100;

[ a, b ] = meshgrid( ( 0 + eps ) : ( 1 - 2 * eps ) / N : ( 1 - eps ), ( 0 + eps ) : ( 0.1 - 2 * eps ) / N : ( 0.1 - eps ) );

[ m, n ] = size( a );

y = zeros( m, n );

z = zeros( m, n );

F0 = [ 0; 0 ];

% load nlsdat1

options=optimset('Display','iter', ...

    'Algorithm', 'trust-region-reflective');

for i = 1 : m 

    for j = 1 : n

        F0 = fsolve( @( x )gzqdyouxia( x, a( i, j ), b( i, j ) ), F0, options );

        y( i, j ) = F0( 1 );

        z( i, j ) = F0( 2 );

    end

end

figure( 'name', '三维图')

subplot( 2, 2, 1 )

surf( a, b, y )

axis equal

title( 'y-a-b' )

subplot( 2, 2, 2 )

surf( a, b, z )

axis equal

title( 'z-a-b' )

% 将a、b、y、z代入方程1，计算误差并画图

diff1 = - y  + z + ( 7999022433351517 / 2305843009213693952 * z + 5509600814111415 / 72057594037927936 * z .^ 3 ) ./ ( 1 - a ) + 1 / 2 * b;

max1 = max( max( diff1 ) )

subplot( 2, 2, 3 )

surf( a, b, diff1 )

title( 'diff1' )

% 将a、b、y、z代入方程2，计算误差并画图

diff2 = - z + y + ( - 4534249478333089 / 576460752303423488 *y - 3611543935826555 / 576460752303423488 * y .^ 3 + 17523 / 50000 * y .^ 5 ) ./ a - 1 / 2 * b;

max2 = max( max( diff2 ) )

subplot( 2, 2, 4 )

surf( a, b, diff2 )

title( 'diff2' )

图：

最大误差：

CODE:

max1 =

      0.000172157751275541

max2 =

     8.02443857179647e-005

17楼2011-01-10 13:02:23

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

xiaowei12

金虫 (正式写手)

应助: 12 (小学生)
金币: 5715
帖子: 530
在线: 217.2小时
虫号: 1675530

精彩的回帖！

18楼2013-10-15 10:13:20

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

crazyiyy

木虫 (小有名气)

应助: 1 (幼儿园)
金币: 1696.5
帖子: 78
在线: 85.6小时
虫号: 2483453

????

19楼2014-07-07 17:18:45

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

相关版块跳转我要订阅楼主 gzqdyouxia 的主题更新