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

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 ];

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' )

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' )

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 ...

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

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

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

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

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

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

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 ...

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 ...

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 ...

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


为什么结果都是零？

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

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

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' )

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  ...

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

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

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


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

关注中。

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

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

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' )

max1 =

      0.000172157751275541


max2 =

     8.02443857179647e-005