MATLAB求解微分方程

求微分方程的解

自牛顿发明微积分以来,微分方程在描述事物运动规律上已发挥了重要的作用。实际应用问题通过数学建模所得到的方程,绝大多数是微分方程。由于实际应用的需要,人们必须求解微分方程。然而能够求得解析解的微分方程十分有限,绝大多数微分方程需要利用数值方法来近似求解。下面介绍如何用 Matlab 来计算微分方程(组)的数值解。

Euler折线法

微分方程的基本数值解法——Euler折线法。

步骤:
分割求解区间,差商代替微商,解代数方程

例子:
在这里插入图片描述解:

等距剖分: a = x 0 < x 1 < x 2 < ⋯ < x n − 1 < x n = b a=x_0<x_1<x_2< \cdots <x_{n-1}<x_n=b a=x0<x1<x2<<xn1<xn=b

步长: h = x k + 1 − x k = ( b − a ) / n , k = 0 , 1 , 2 , . . . , n − 1 h=x_{k+1}-x_k=(b-a)/n,k=0,1,2,...,n-1 h=xk+1xk=(ba)/n,k=0,1,2,...,n1

y ( x k + 1 ) ≈ y ( x k ) + h y ′ ( x k ) y(x_{k+1})\approx y(x_k)+hy\prime(x_k) y(xk+1)y(xk)+hy(xk)

取步长 h = (2 - 0)/n = 2/n,得差分方程
在这里插入图片描述
当 h=0.4,即 n=5 时,代码如下:

clear
a=0;  b=2;
h=0.4;
n=(b-a)/h+1; 
x(1)=0;  y(1)=1;
for i=1:n-1        
    y(i+1)=y(i)+h*(y(i)+2*x(i)/y(i)^2);
    x(i+1)=x(i)+h;
end
y1=((5/3)*exp(3*x)-2*x-2/3).^(1/3);
plot(x,y1,x,y);
h = legend('解析解','Euler解');

在这里插入图片描述

解得:
y = ( 5 3 e 3 x − 2 x − 2 3 ) 1 / 3 y=(\frac{5}{3}e^{3x}-2x-\frac{2}{3})^{1/3} y=(35e3x2x32)1/3

Runge-Kutta方法

龙格-库塔法(Runge-Kutta methods)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。

同样是上面的例子:
在这里插入图片描述
在这里插入图片描述

clear;
a=0; b=2; h=0.4;
n=(b-a)/h+1; 
x(1)=0; y(1)=1;  
for i=1:n-1
    L1=y(i)+2*x(i)/y(i)^2;
    t=x(i)+h/2;t1=y(i)+h*L1/2;
    L2=t1+2*t/t1^2;t1=y(i)+h*L2/2;
    L3=t1+2*t/t1^2;t1=y(i)+h*L3/2;
    L4=t1+2*(x(i)+h)/t1^2;
    y(i+1)=y(i)+h*(L1+2*L2+2*L3+L4)/6; 
    x(i+1)=x(i)+h;
end
plot(x,y)

在这里插入图片描述

MATLAB自带的ODE求解器

ODE (常微分方程)

没有一种算法可以有效地解决所有的 ODE 问题,因此MATLAB 提供了多种ODE求解器,对于不同的ODE,可以调用不同的求解器。
在这里插入图片描述

基本语法
[T,Y] = solver(odefun,tspan,y0)

其中 y0 为初值条件,tspan为求解区间;Matlab在数值求解时自动对求解区间进行分割,T (向量) 中返回的是分割点的值(自变量),Y (向量) 中返回的是解函数在这些分割点上的函数值。
solver 为Matlab的ODE求解器(可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb)

posted @   Patrick-Rex  阅读(59)  评论(0编辑  收藏  举报  
相关博文:
阅读排行:
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 全网最简单!3分钟用满血DeepSeek R1开发一款AI智能客服,零代码轻松接入微信、公众号、小程
· .NET 10 首个预览版发布,跨平台开发与性能全面提升
点击右上角即可分享
微信分享提示