数值分析 常微分方程数值解 欧拉显式方法 预估-校正法 4阶龙格-库塔法

一、题目要求

第3次作业：用有限差分法求解下面的定解问题

要求：利用欧拉显式方法，预估-校正法，4阶龙格-库塔法；分析h取不同的值对计算结果的影响。

二、原理分析

2.1 欧拉显式方法

2.2预估-校正法

      

2.3 4阶龙格-库塔法

三、代码实现

 代码用java语言实现

因为折线图搞得有点头皮发麻，代码写的比较乱

import org.jfree.chart.ChartFactory;
import org.jfree.chart.ChartFrame;
import org.jfree.chart.JFreeChart;
import org.jfree.chart.plot.PlotOrientation;
import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;

import javax.swing.*;

public class Main {
    public static void main(String[] args) {
        /**
         * 数值分析
         */
        Run run = new Run();
        run.programImplementation();
    }
}
class Run{
    //创建一个数据集
    XYSeriesCollection dataset = new XYSeriesCollection();
    XYSeriesCollection dataset1 = new XYSeriesCollection();
    XYSeriesCollection dataset2 = new XYSeriesCollection();

    public void programImplementation() {
        //设计基本参数
        double x0 = 0.0;
        double y0 = 1.0;
        double h1 = 0.5;
        double h2 = 0.2;
        double h3 = 0.1;
        double h4 = 0.05;



        explicitEuler(x0, y0, h1);//创建四条数据链，迭代计算并赋值
        explicitEuler(x0, y0, h2);
        explicitEuler(x0, y0, h3);
        explicitEuler(x0, y0, h4);



        // 利用数据集创建一个折线图
        JFreeChart chart = ChartFactory.createXYLineChart(
                "explicitEuler", // 图表标题
                "X", // X轴标签
                "Y", // Y轴标签
                dataset, // 数据集
                PlotOrientation.VERTICAL, // 图表方向
                true, // 是否显示图例
                true, // 是否生成工具提示
                false // 是否生成URL链接
        );




        // 创建一个图形窗口
        ChartFrame frame = new ChartFrame("折线图示例", chart);
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame.pack();
        frame.setVisible(true);


        predictorCorrectorMethod(x0, y0, h1);//创建四条数据链，迭代计算并赋值

        predictorCorrectorMethod(x0, y0, h2);

        predictorCorrectorMethod(x0, y0, h3);

        predictorCorrectorMethod(x0, y0, h4);

        JFreeChart chart1 = ChartFactory.createXYLineChart(
                "predictorCorrectorMethod", // 图表标题
                "X", // X轴标签
                "Y", // Y轴标签
                dataset1, // 数据集
                PlotOrientation.VERTICAL, // 图表方向
                true, // 是否显示图例
                true, // 是否生成工具提示
                false // 是否生成URL链接
        );
        // 创建一个图形窗口
        ChartFrame frame1 = new ChartFrame("折线图示例1", chart1);
        frame1.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame1.pack();
        frame1.setVisible(true);

        fourThorderRungeKuttaMethod(x0, y0, h1);//创建四条数据链，迭代计算并赋值

        fourThorderRungeKuttaMethod(x0, y0, h2);

        fourThorderRungeKuttaMethod(x0, y0, h3);

        fourThorderRungeKuttaMethod(x0, y0, h4);

        JFreeChart chart2 = ChartFactory.createXYLineChart(
                "fourThorderRungeKuttaMethod", // 图表标题
                "X", // X轴标签
                "Y", // Y轴标签
                dataset2, // 数据集
                PlotOrientation.VERTICAL, // 图表方向
                true, // 是否显示图例
                true, // 是否生成工具提示
                false // 是否生成URL链接
        );
        // 创建一个图形窗口
        ChartFrame frame2 = new ChartFrame("折线图示例2", chart2);
        frame2.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        frame2.pack();
        frame2.setVisible(true);

    }


    public void fourThorderRungeKuttaMethod(double x, double y, double h) {
        XYSeries series;
        if (h == 0.5) {
            series = new XYSeries("fourThorderRungeKuttaMethod(h=0.5)");
        } else if (h == 0.2) {
            series = new XYSeries("fourThorderRungeKuttaMethod(h=0.2)");
        } else if (h == 0.1) {
            series = new XYSeries("fourThorderRungeKuttaMethod(h=0.1)");
        }else {
            series = new XYSeries("fourThorderRungeKuttaMethod(h=0.05)");
        }
        int s = (int) (1.0 / h);
        for (int i = 0; i < s; i++) {//计算4阶龙格-库塔法
            series.add(x,y);
            double k1 = Math.exp(-(x + y)) * Math.sin(x * y);
            double xk = x + 0.5 * h;
            double yk = y + 0.5 * h * k1;
            double k2 = Math.exp(-(xk + yk)) * Math.sin(xk * yk);
            yk = y + 0.5 * h * k2;
            double k3 = Math.exp(-(xk + yk)) * Math.sin(xk * yk);
            xk = x + h;
            yk = y + h * k3;
            double k4 = Math.exp(-(xk + yk)) * Math.sin(xk * yk);
            y = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
            x = x + h;

            System.out.print("4阶龙格-库塔方法(h =" + h + ")：第" + (i + 1) + "次迭代结果为:" + y + "\t");
        }
        series.add(x,y);
        System.out.println();
        //将数据链数据输入数据集
        dataset2.addSeries(series);
    }

    public void predictorCorrectorMethod(double x, double y, double h) {//计算预估-校正法
        XYSeries series;
        if (h == 0.5) {
            series = new XYSeries("predictorCorrectorMethod(h=0.5)");
        } else if (h == 0.2) {
            series = new XYSeries("predictorCorrectorMethod(h=0.2)");
        } else if (h == 0.1) {
            series = new XYSeries("predictorCorrectorMethod(h=0.1)");
        }else {
            series = new XYSeries("predictorCorrectorMethod(h=0.05)");
        }
        int s = (int) (1.0 / h);
        for (int i = 0; i < s; i++) {
            series.add(x,y);
            double k1 = Math.exp(-(x + y)) * Math.sin(x * y);
            double xk = x + h;
            double yk = y + h * k1;
            double k2 = Math.exp(-(xk + yk)) * Math.sin(xk * yk);
            y = y + (h / 2) * (k1 + k2);
            x = x + h;

            System.out.print("预估-校正法(h =" + h + ")：第" + (i + 1) + "次迭代结果为:" + y + "\t");
        }
        series.add(x,y);
        System.out.println();
        //将数据链数据输入数据集
        dataset1.addSeries(series);
    }


    public void explicitEuler(double x, double y, double h) {//计算欧拉显式方法
        XYSeries series;
        if (h == 0.5) {
            series = new XYSeries("explicitEuler(h=0.5)");
        } else if (h == 0.2) {
            series = new XYSeries("explicitEuler(h=0.2)");
        } else if (h == 0.1) {
            series = new XYSeries("explicitEuler(h=0.1)");
        }else {
            series = new XYSeries("explicitEuler(h=0.05)");
        }
        int s = (int) (1.0 / h);

        for (int i = 0; i < s; i++) {
            series.add(x,y);
            y = y + h * (Math.exp(-(x + y)) * Math.sin(x * y));
            x = x + h;

            System.out.print("显式欧拉方法(h =" + h + ")：第" + (i + 1) + "次迭代结果为:" + y + "\t");
        }
        series.add(x,y);
        System.out.println();
        //将数据链数据输入数据集
        dataset.addSeries(series);
    }
}

四、结果展示

根据步长的不同，拥有不同的迭代次数

当h = 0.5时，两次迭代

当 h = 0.2时，五次迭代

当 h = 0.1时，十次迭代

当 h = 0.05时，二十次迭代

4.1欧拉显式方法

欧拉显式方法，不同步长迭代结果如下：

4.2预估-校正法

预估-校正法，不同步长迭代结果如下：

4.3 4阶龙格-库塔法

4阶龙格-库塔法，不同步长迭代结果如下：

4.4 结果分析

通过上述结果，可以看出4阶龙格-库塔法的效果最好，其对步长的鲁棒性更好，步长大小对于结果的影响不大，欧拉显式方法效果最差，不同的步长对结果影响很大。