求微分方程组dsolve('Dx=-x+y^2','Dy=-2*y+x^2','t')的详细代码,还要画出xy的函数图像

求微分方程组dsolve('Dx=-x+y^2','Dy=-2*y+x^2','t')的详细代码,还要画出xy的函数图像
求微分方程组dsolve('Dx=-x+y^2','Dy=-2*y+x^2','t')的详细代码,还要画出xy的函数图像

求微分方程组dsolve('Dx=-x+y^2','Dy=-2*y+x^2','t')的详细代码,还要画出xy的函数图像
!
!pm_test.f90
!lcx-fortran
!
!Created by apple on 12-2-3.
!Copyright 2012 apple. All rights reserved.
!
module rkf_45
implicit none
contains
subroutine F(t,n,x,dx)
implicit none
integer n
real*8 t
real*8,dimension(n) :: x,dx
dx(1) = x(1)+x(2)**2
dx(2) = x(1)**2+2*x(2)
return
end subroutine
subroutine rkf4(t,a,b,xa,n,M,z)
implicit none
real*8 a,b
real*8 h
real*8,dimension(n) :: k1,k2,k3,k4
integer n,m
integer j
real*8,dimension(n) :: xa,dz
real*8,dimension(m+1) :: t
real*8,dimension(M+1,n) :: z
h=(b-a)/m
t(1)=a
z(1,:)=xa
do j=1,m
call f(t(j),n,z(j,:),dz)
k1=dz*h
call f(t(j)+h/2.0d0,n,z(j,:)+k1/2.0d0,dz)
k2=dz*h
call f(t(j)+h/2.0d0,n,z(j,:)+k2/2.0d0,dz)
k3=dz*h
call f(t(j)+h,n,z(j,:)+k3,dz)
k4=dz*h
z(j+1,:)=z(j,:)+(k1+2*k2+2*k3+k4)/6.0d0
t(J+1)=t(j)+h
end do
return
end subroutine
end module
program main
use rkf_45
implicit none
integer,parameter :: n=2
real*8,dimension(n) :: xa
real*8 a,b
integer m
real*8,allocatable :: t(:)
real*8,allocatable :: z(:,:)
integer i
integer,parameter :: fileid = 10
character(len=20) :: filename = "rk45.txt"
open(fileid,file=filename)
a=0.0d0 !初始时刻
b=0.2d0 !结束时刻
m=100d0 !步数
allocate(t(m+1))
allocate(z(M+1,n))
xa=(/6.0,4.0/) ! 初值
call rkf4(t,a,b,xa,n,M,z)
do i=1,m+1
write(fileid,"(2(2x,f17.14))") z(i,:)
end do
stop
end
至于图像,将生成的txt文件导入作图软件就可以显示x—y图像了,由于你没给任何初值,积分区间,我只能事先拟定了一组参数.这是fortran语言编写,不知你是用c,pascal,还是其他语言.这是数值法rk4求非线性微分方程,对于线性的,不用编程就可以做.