隐式差分格式怎么求解

cuixindong

言深深 发表于 2013-3-28 18:16

                               
登录/注册后可看大图

额，汗···帖子已经好些日子了吧···

我现在在模拟一个方程组，包含声波和重力波，但是解一直增长的，计算不稳定。我试过很多办法，都没用，想找一个程序参考一下。不知楼主可否提一下建议？

言深深

cuixindong 发表于 2013-3-28 19:56

                               
登录/注册后可看大图

我现在在模拟一个方程组，包含声波和重力波，但是解一直增长的，计算不稳定。我试过很多办法，都没用，想 ...

数值方法求解的时候就是容易出现这样的情况，建议你看看程序有没有编错，另外就是可以试试一些收敛性比较好的办法，比如四阶龙格库塔法，误差比较小。不放查看依稀数值方法的书，上面有些例子是可以借鉴的···祝你成功

cuixindong

言深深 发表于 2013-3-28 20:46

                               
登录/注册后可看大图

数值方法求解的时候就是容易出现这样的情况，建议你看看程序有没有编错，另外就是可以试试一些收敛性比较 ...

版主，谢谢一直回复和帮助。您刚刚提到的这个依稀数值方法是什么？

言深深

cuixindong 发表于 2013-3-29 07:55

                               
登录/注册后可看大图

版主，谢谢一直回复和帮助。您刚刚提到的这个依稀数值方法是什么？

额，汗，打错了··是一些···

cuixindong

版主，您有关于四阶龙格库塔法的材料吗？能发一下吗？

cuixindong

言深深 发表于 2013-3-29 08:39

                               
登录/注册后可看大图

额，汗，打错了··是一些···

版主，您有关于四阶龙格库塔法的材料吗？能发一下吗？

言深深

cuixindong 发表于 2013-3-29 09:39

                               
登录/注册后可看大图

版主，您有关于四阶龙格库塔法的材料吗？能发一下吗？

我去···帅哥，程序都是针对各人的问题的···
这是朋友写针对ta的薛定谔方程写的，你肯定用不了···但是可以让你了解一下这个方法啥意思···
do g=1,mm
  sum=0
  aaa=0
  do i=1,N
          s(g,i)=u1(i)*conjg(u1(i+1))-conjg(u1(i))*u1(i+1)
        sum=sum+s(g,i)
        aaa=aaa+abs(u1(i))
  end do

  do i=1,N
    if (i==ll) then
        u2(i)=h/4.0*c*((u1(i-1)+u1(i+1))/2.0-(e+gm*(u1(i))**r*(conjg(u1(i)))**r)*u1(i))+u1(i)
        else
    u2(i)=h/4.0*c*((u1(i-1)+u1(i+1))/2.0-(gm*(u1(i))**r*(conjg(u1(i)))**r)*u1(i))+u1(i)
        end if
  end do
    u2(0)=u2(N)
        u2(N+1)=u2(1)
  
  do i=1,N
    if (i==ll) then
        u3(i)=h/3.0*c*((u2(i-1)+u2(i+1))/2.0-(e+gm*(u2(i))**r*(conjg(u2(i)))**r)*u2(i))+u1(i)
        else
    u3(i)=h/3.0*c*((u2(i-1)+u2(i+1))/2.0-(gm*(u2(i))**r*(conjg(u2(i)))**r)*u2(i))+u1(i)
   
        end if
  end do
   u3(0)=u3(N)
   u3(N+1)=u3(1)

  do i=1,N
    if (i==ll) then     
        u4(i)=h/2.0*c*((u3(i-1)+u3(i+1))/2.0-(e+gm*(u3(i))**r*(conjg(u3(i)))**r)*u3(i))+u1(i)
        else
    u4(i)=h/2.0*c*((u3(i-1)+u3(i+1))/2.0-(gm*(u3(i))**r*(conjg(u3(i)))**r)*u3(i))+u1(i)
   
        end if
  end do
     u4(0)=u4(N)
         u4(N+1)=u4(1)

  do i=1,N
    if (i==ll) then
        u5(i)=h*c*((u4(i-1)+u4(i+1))/2.0-(e+gm*(u4(i))**r*(conjg(u4(i)))**r)*u4(i))+u1(i)
        else
        u5(i)=h*c*((u4(i-1)+u4(i+1))/2.0-(gm*(u4(i))**r*(conjg(u4(i)))**r)*u4(i))+u1(i)  
   
        end if
  end do
     u5(0)=u5(N)
         u5(N+1)=u5(1)
  
  do  i=0,N+1   
    u1(i)=u5(i)
  end do
enddo

ps：改程序不针对隐士差分格式，仅回答问题····

cuixindong

言深深 发表于 2013-3-29 09:47

                               
登录/注册后可看大图

我去···帅哥，程序都是针对各人的问题的···
这是朋友写针对ta的薛定谔方程写的，你肯定用不了·· ...

谢谢版主。

eeecho

苏叶绿阳

来此向大家学习了