На главную Наши проекты:
Журнал   ·   Discuz!ML   ·   Wiki   ·   DRKB   ·   Помощь проекту
ПРАВИЛА FAQ Помощь Участники Календарь Избранное RSS
msm.ru
! Друзья, соблюдайте, пожалуйста, правила форума и данного раздела:
Данный раздел не предназначен для вопросов и обсуждений, он содержит FAQ-заготовки для разных языков программирования. Любой желающий может разместить здесь свою статью. Вопросы же задавайте в тематических разделах!
• Если ваша статья может быть перенесена в FAQ соответствующего раздела, при условии, что она будет оформлена в соответствии с Требованиями к оформлению статей.
• Чтобы остальным было проще понять, указывайте в описании темы (подзаголовке) название языка в [квадратных скобках]!
Модераторы: Модераторы
  
> Метод Милна для системы ОДУ , [Pascal]
    ExpandedWrap disabled
      Program Miln;
      Uses Crt;
      var
        x0,y0,z0,h:real;
        k1,k2,k3,k4:real;
        r1,r2,r3,r4:real;
        eps,abs_pogr:real;
       
        z,zkor,zpr,ypr,ykor,x,y:array [0..10] of real;
        i:integer;
       
      function f1(xa,ya,za:real):real;
      begin
        f1:=2*sqr(xa)+2*ya+za;
      end;
       
      function f2(xa,ya,za:real):real;
      begin
        f2:=1-2*sqr(xa)+2*ya-za;
      end;
       
      begin
        Clrscr;
        eps:=10e-6;
        x0:=0;
        y0:=1;
        z0:=1;
        h:=0.1;
        x[0]:=x0;
        y[0]:=y0;
        z[0]:=y0;
        i:=0;
        { Gotovim 1-e 3 tochki po metodu runge-kutta }
        while i<=3  do
          begin
            k1:=h*f1(x[i],y[i],z[i]);
            r1:=h*f2(x[i],y[i],z[i]);
            k2:=h*f1(x[i]+h/2,y[i]+k1/2,z[i]+r1/2);
            r2:=h*f2(x[i]+h/2,y[i]+k1/2,z[i]+r1/2);
            k3:=h*f1(x[i]+h/2,y[i]+k2/2,z[i]+r2/2);
            r3:=h*f2(x[i]+h/2,y[i]+k2/2,z[i]+r2/2);
            k4:=h*f1(x[i]+h,y[i]+k3,z[i]+r3);
            r4:=h*f2(x[i]+h,y[i]+k3,z[i]+r3);
       
            y[i+1]:=y[i]+(k1+2*k2+2*k3+k4)/6;
            z[i+1]:=z[i]+(r1+2*r2+2*r3+r4)/6;
       
            x[i+1]:=x[i]+h;
            i:=i+1;
          end;
       
        i:=4;
        while x[i]<=1.0+h{eps> abs_pogr} do
          begin
        { etap prognoza i korrektsii}
            ypr[i]:=y[i-4]+(4*h)/3*(2*f1(x[i-3],y[i-3],z[i-3])-f1(x[i-2],y[i-2],z[i-2])+2*f1(x[i-1],y[i-1],z[i-1]));
       
            ykor[i]:=y[i-2]+(h/3)*(f1(x[i-2],y[i-2],z[i-2])+4*f1(x[i-1],y[i-1],z[i-1])+f1(x[i],ypr[i],z[i]));
       
            zpr[i]:=z[i-4]+(4*h)/3*(2*f2(x[i-3],y[i-3],z[i-3])-f2(x[i-2],y[i-2],z[i-2])+2*f2(x[i-1],y[i-1],z[i-1]));
       
            zkor[i]:=z[i-2]+(h/3)*(f2(x[i-2],y[i-2],z[i-2])+4*f2(x[i-1],y[i-1],z[i-1])+f2(x[i],zpr[i],z[i]));
       
       
            abs_pogr:=abs(ykor[i]-ypr[i])/29;
            if abs_pogr>eps then
              y[i]:=ykor[i]
            else
                y[i]:=ypr[i];
            abs_pogr:=abs(zkor[i]-zpr[i])/29;
            if abs_pogr>eps then
              z[i]:=zkor[i]
            else
                z[i]:=zpr[i];
            x[i+1]:=x[i]+h;
            i:=i+1;
          end;
        WriteLn;
        for i:=0 to 10 do
          begin
            WriteLn(x[i]:10:4,' ',y[i]:10:4,' ',z[i]:10:4);
          end;
        Readln;
      end.
    Сообщение отредактировано: Jin X -
      Было бы совсем неплохо словами описать алгоритм, так как он не очевиден... Для меня, по крайней мере...
        http://window.edu.ru/window_catalog/files/r18677/Mtduksi5.pdf
        0 пользователей читают эту тему (0 гостей и 0 скрытых пользователей)
        0 пользователей:


        Рейтинг@Mail.ru
        [ Script execution time: 0,0194 ]   [ 16 queries used ]   [ Generated: 29.03.24, 15:57 GMT ]