[½»Á÷] ¡¾ÇóÖú¡¿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:

²»ºÃÒâ˼£¬Î´ÄÜÉÆʼÉÆÖÕ£¬ºÇºÇ¡£

²»ÖªµÀÕâÊDz»ÊÇÕýÈ·µÄ£º

[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