• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    迪恩网络公众号

数值分析实验之常微分方程的数值解法(MATLAB实现)

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

一、实验目的

科学技术中常常要求解常微分方程的定解问题,所谓数值解法就是求未知函数在一系列离散点处的近似值。

二、实验原理

三、实验程序

  1. 尤拉公式程序

     

 

      2、3、4的尤拉公式的程序参上改写。

 四、实验内容

   

 

五、实验代码及运行结果

    • MATLAB代码:

定义函数:
function [A1,A2,B1,B2,C1,C2]=euler23(a,b,n,y0)
%欧拉法解一阶常微分方程
%初始条件y0
h = (b-a)/n; %步长h
%区域的左边界a
%区域的右边界b
x = a:h:b;
m=length(x);
%前向欧拉法
y = y0;
for i=2:m
    y(i)=y(i-1)+h*oula(x(i-1),y(i-1));
    A1(i)=x(i);
    A2(i)=y(i);
end
plot(x,y,\'r-\');
hold on;
%改进欧拉法
y = y0;
for i=2:m
    y(i)= y(i-1)+h/2*( oula(x(i-1),y(i-1))+oula(x(i),y(i-1))+h*(oula(x(i-1),x(i-1))));
    B1(i)=x(i);
    B2(i)=y(i);
end
plot(x,y,\'m-\');
hold on;
%欧拉两步公式
y=y0;
y(2)=y(1)+h*oula(x(1),y(1));
for i=2:m-1
    y(i+1)=y(i-1)+2*h*oula(x(i),y(i));
    C1(i)=x(i);
    C2(i)=y(i);
end
plot(x,y,\'b-\');
hold on;
%精确解用作图
xx = x;
f = dsolve(\'Dy=-y+1\',\'y(0)=0\',\'x\');%求出解析解
y = subs(f,xx); %将xx代入解析解,得到解析解对应的数值
plot(xx,y,\'k--\');
legend(\'前向欧拉法\',\'改进欧拉法\',\'欧拉两步法\',\'解析解\');

function f=oula(x,y)
f=-y+1;

命令行窗口:
[A1,A2,B1,B2,C1,C2]=euler23(0,1,10,0)

  运行结果:

      

  N=50时:

      

 

N=100时:

      

 

   故得精度越大时,几种方法求解值与准确值越来越接近。

 

     • 另解

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 0;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
xx = x;
f = dsolve(\'Dy=-y+1\',\'y(0)=0\',\'x\');
z  = subs(f,xx); 
y(2,:) = z;
plot(x, y);



function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = 1-y0;
        d = y0 + h*f0;
        y1 = calculate(y0, x1, d, h);
        %result = calculate(x1, d, h);
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end

function result = calculate(y0, x1, y1, h)

    acc = -6;
    now = 0.0;
    z1 = y1;

    while now >= -6

        z0 = z1;
        f0 =1-z0;
        z1 = y0 + h*f0;
        now = log10(abs(z1-z0));

    end
    result = z1;
end

  运行结果:

     

 


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
DELPHITDownLoadURL下载网络文件发布时间:2022-07-18
下一篇:
Delphi自带的字符串分割函数split发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap