Fortran | 普通差分法解非线性常微分方程组（洛伦兹对流方程组）

本文链接：https://blog.youkuaiyun.com/qq_53185200/article/details/127561592

0.写在前面

今天上机，老师直接让解方程组，我们上节课只教了龙格-库塔解法，我直接一脸懵逼......好在想到了，就是用最一般的差分解。

1.待求方程组

$\\\frac{\mathrm{d} x }{\mathrm{d} t} =\sigma (y-x) \\ \frac{\mathrm{d} y }{\mathrm{d} t} =x(\rho -z)-y\\ \frac{\mathrm{d} z }{\mathrm{d} t} =xy-\beta z$

2.初始条件和要求

σ=10，ρ=28，β=8.0/3

t0=0，x0=0，y0=1，z0=0

Δt=0.01 N=3000

3.代码

program exp1
    implicit none
    real x0,x(1:3000),y0,y(1:3000),z0,z(1:3000),t0,t(1:3000),ts,a,b,c
    integer N,i
    ts=0.01
    n=3000
    x0=0.0
    y0=1.0
    z0=0.0
    t0=0.001
    a=10.0
    b=28.0
    c=8.0/3
    x(1)=x0+ts*a*(y0-x0)
    y(1)=y0+ts*x0*(b-z0)-y0*ts
    z(1)=z0+ts*(x0*y0-c*z0)
    t(1)=t0+ts
    do i=2,n
        x(i)=x(i-1)+ts*a*(y(i-1)-x(i-1))
        y(i)=y(i-1)+ts*x(i-1)*(b-z(i-1))-y(i-1)*ts 
        z(i)=z(i-1)+ts*(x(i-1)*y(i-1)-c*z(i-1))
        t(i)=t(i-1)+ts
    enddo
    open(1,file='./output_exp1.csv')
    write(1,'(3(f13.8,a),f13.8)') t0,',',x0,',',y0,',',z0
    do i=1,n
        write(1,'(3(f13.8,a),f13.8)') t(i),',',x(i),',',y(i),',',z(i)
    enddo
    close(1)
end program exp1

#因为可视化的代码是老师给的，不好放出来，大家看个效果就行

#exp2的初值z0改为0.001，其他条件不变