Потск корней методом Мюллера

Форум пользователей пакета Maple

Модератор: Admin

Henry_12
Сообщения: 4
Зарегистрирован: Ср дек 14, 2005 6:22 pm
Контактная информация:

Потск корней методом Мюллера

Сообщение Henry_12 » Ср дек 14, 2005 6:34 pm

Я уже запарился, процедура каждый раз выдает разнее ответы и иногда редко находит правильный ответ, в чем проблемма мож кто нить что нить увидит чего я не вижу

Код: Выделить всё

restart;
> RootM:=proc(F,eps,x)
>          local k1,NewZ,s0,s1,s2;
>          print("******* Начальное приближение ***********");
>          s0:=x[1]+I*x[2];
>          print(s0);
>          s1:=s0-Complex(0.1,0.1);
>          print(s1);
>          s2:=s0+Complex(0.1,0.1);
>          print(s2);
>          print("*************** Ответ *******************");
>          k1:=0;
>          NewZ:=CalcD(s0,s1,s2,F);
>         
> while (abs(NewZ-s2)>eps)
>            do
> #print(abs(NewZ-s2));
>              k1:=k1+1;
>              s0:=s1;
>              s1:=s2;
>              s2:=NewZ;
>              NewZ:=CalcD(s0,s1,s2,F);
>            end do;
>          print(NewZ);
>          #print(k);
>        end proc:
> ###############################################################################################
> CalcD:=proc(Z0,Z1,Z2,fun)
>          local zc1,zc2,g1,g2,g3,k,D,Distance1,Distance2,zc,c1,c2,c3; 
>          #print(Z0,Z1,Z2);
>          g1:=c1*Z0^2+c2*Z0+c3=fun(Z0);
>          g2:=c1*Z1^2+c2*Z1+c3=fun(Z1);
>          g3:=c1*Z2^2+c2*Z2+c3=fun(Z2);
>            #print(fun(Z0));
>            #print(g1);
>            #print(g2);
>            #print(g3);
>          k:=fsolve({g1,g2,g3},{c1,c2,c3},complex);
>            #print(k[1]);
>            #print(rhs(k[2]));
>            #print(k[3]);
>            #print(k[2]);
>          c1:=rhs(k[1]); c2:=rhs(k[2]); c3:=rhs(k[3]);
>          D:=c1^2-4*c2*c3;
>          #print(D);
>          zc1:=(-rhs(k[1])+sqrt(D))/2;
>          zc2:=(-rhs(k[1])-sqrt(D))/2;
>          #print(rhs(k[2]));
>          #print(zc2);
>          Distance1:=Z2-zc1;
>          Distance2:=Z2-zc2;
>          zc:=zc1;
>          if (Distance1>Distance2) then zc:=zc2 end if;
>          print(zc);
>          return(zc);
>        end proc:
>f:=t->t^2+1;

                                      2
                           f := t -> t  + 1


> s:=array([1.0,1.0]);

                           s := [1.0, 1.0]

> RootM(f,0.001,s);

Кто нить проверьте у себя и посоветуйте что нибудь пожалуйста. Или если у кого есть подобная процедурка вычисления корней методом Мюллера буду благодарен:)

Break
Сообщения: 159
Зарегистрирован: Вс окт 09, 2005 2:10 am
Откуда: Петербург

Сообщение Break » Чт дек 15, 2005 6:33 pm

Возможно, в ваших иследованиях, будет полезна эта ссылка:
http://exponenta.krasu.ru/educat/system ... _1.asp.htm

Но там Маткад.

Henry_12
Сообщения: 4
Зарегистрирован: Ср дек 14, 2005 6:22 pm
Контактная информация:

Сообщение Henry_12 » Чт дек 15, 2005 8:16 pm

Спосибо конечно, но сам метод вплане математики меня не интересует, информация по нему есть. Меня интересует именно реализация и именно в пакете Maple

Break
Сообщения: 159
Зарегистрирован: Вс окт 09, 2005 2:10 am
Откуда: Петербург

Сообщение Break » Чт дек 15, 2005 8:22 pm

Там еще были размышления на тему возникновения ложных корней в зависимости от точности приближения... Просто исправлять чужие алгоритмы - дело не благодарное! :)

Break
Сообщения: 159
Зарегистрирован: Вс окт 09, 2005 2:10 am
Откуда: Петербург

Сообщение Break » Чт дек 15, 2005 9:14 pm

Пока нашел одну ошибку:

Вы пишите,

c1:=rhs(k[1]); c2:=rhs(k[2]); c3:=rhs(k[3]);

- так не надо, т.к. каждый раз Мапл сортирует решения в произвольном порядке. Убедитесь, если напишите print("k",k): после fsolve.

Пишите assign(k). Можете проверить - он работает корректно:

print("k",k):
assign(k);
print(c1,c2,c3) ;

Henry_12
Сообщения: 4
Зарегистрирован: Ср дек 14, 2005 6:22 pm
Контактная информация:

Сообщение Henry_12 » Пт дек 16, 2005 4:24 pm

Проверил, так не рабоает: на строчке print(c1,c2,c3) вылетает Error, (in CalcD) cannot evaluate boolean:

Break
Сообщения: 159
Зарегистрирован: Вс окт 09, 2005 2:10 am
Откуда: Петербург

Сообщение Break » Пт дек 16, 2005 4:50 pm

Какая версия Мапла?

Но, хоть убедились, что k каждый раз в произвлольном порядке выдает? Это факт. Только ?assign помогает.
Приведу, на всяк случ. свой вариант:

Henry_12
Сообщения: 4
Зарегистрирован: Ср дек 14, 2005 6:22 pm
Контактная информация:

Сообщение Henry_12 » Пт дек 16, 2005 7:07 pm

Maple 6 :)
C assign понятно. Но дело тут по большому счету не в нем