24小时热门版块排行榜    

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

gzqdyouxia

金虫 (著名写手)


[交流] 【求助】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 ]
回复此楼

» 猜你喜欢

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

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

查看全部散金贴

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

blue.wh

铁杆木虫 (知名作家)


引用回帖:
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的回帖
查看全部 19 个回答
★ ★ ★
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 ]
2楼2010-12-27 19:40:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

gzqdyouxia

金虫 (著名写手)


引用回帖:
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的回帖
★ ★
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)
就可以了啊
4楼2010-12-27 22:56:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +3 DpSrDtM079iu 2026-08-22 3/150 2026-08-23 19:55 by lXgU0RvNTDC9
[找工作] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 2JOx3r2CYEgw 2026-08-22 5/250 2026-08-23 19:43 by lXgU0RvNTDC9
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 2JOx3r2CYEgw 2026-08-22 4/200 2026-08-23 19:34 by lXgU0RvNTDC9
[博后之家] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 2JOx3r2CYEgw 2026-08-22 6/300 2026-08-23 19:31 by lXgU0RvNTDC9
[论文投稿] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +3 2JOx3r2CYEgw 2026-08-21 8/400 2026-08-23 18:58 by lXgU0RvNTDC9
[找工作] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +3 h4CP7TrQR8Lg 2026-08-22 3/150 2026-08-23 16:30 by lXgU0RvNTDC9
[论文投稿] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +3 DpSrDtM079iu 2026-08-22 3/150 2026-08-23 15:33 by KM5EcsNQRBPn
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 2JOx3r2CYEgw 2026-08-21 5/250 2026-08-23 14:33 by KM5EcsNQRBPn
[基金申请] 让我中一个面上吧! +12 大萍1987 2026-08-20 14/700 2026-08-23 10:30 by wrm
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +3 2JOx3r2CYEgw 2026-08-22 5/250 2026-08-23 07:41 by OEbVnUOu01ol
[考博] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 2JOx3r2CYEgw 2026-08-21 10/500 2026-08-23 07:17 by OEbVnUOu01ol
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 2JOx3r2CYEgw 2026-08-22 7/350 2026-08-23 03:55 by OEbVnUOu01ol
[论文投稿] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 2JOx3r2CYEgw 2026-08-22 5/250 2026-08-23 03:52 by OEbVnUOu01ol
[论文投稿] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 QTy3jDtz1uLt 2026-08-21 9/450 2026-08-23 03:04 by OEbVnUOu01ol
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 QTy3jDtz1uLt 2026-08-21 11/550 2026-08-23 02:55 by OEbVnUOu01ol
[基金申请] 今天基金会出结果吗?20260819 +16 kkkl_v 2026-08-19 17/850 2026-08-22 16:12 by 阿布Abu
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 应该是下周三26日公布了吧? +4 哈哈蛤? 2026-08-21 4/200 2026-08-21 10:58 by Vivilian
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
信息提示
请填处理意见