请问怎么用matlab求这个微分方程的数值解并绘制解函数的图形? |
您所在的位置：网站首页 › 中国华融股份董事长是谁 › 请问怎么用matlab求这个微分方程的数值解并绘制解函数的图形? |
请问怎么用matlab求这个微分方程的数值解并绘制解函数的图形?
|
摘要：讲解欧拉法求解微分方程原理，通过MATLAB程序求解实例。 求解微分方程的时候，如果不能将求出结果的表达式，则可以对利用数值积分对微分方程求解，获取数值解。欧拉方法是最简单的一种数值解法。本文理论部分来自知乎作者云端之下的文章“常微分方程——数值解——欧拉方法”，文章链接为： https://zhuanlan.zhihu.com/p/464118275  实例求解微分方程dy/dt=-y+t+1,y(0)=1,t的取值为0到2，步长h=0.1,用欧拉法求解微分方程并将结果与y(t)=exp(-t)+t比较。 主程序 clc; clear all; close all; h = 0.1;%步长 y0 = 1;%初值 t = 0:h:2;%x范围 y = exp(-t)+t;%真解 n = length(t); numy = zeros(1,n); f2 = -y0; numy(1) = y0; %欧拉法计算 for i=2:n numy(i) = euler1(y0,h,f2); y0 = numy(i); f2 = f1(t(i),y0); end %绘图 figure; plot(t,y,'r-','linewidth',1); hold on; plot(t,numy,'b-','linewidth',1); xlabel('t'); grid on; title('Euler法求系统的输出响应'); ylabel('输出响应y(t)'); legend('真解','Euler法','location','northwest'); wucha_euler = (numy-y).^2; disp('Euler法误差平方:'); wucha_euler自定义函数euler1.m function y = euler1(y0,h,f2) %%输入参数 y0表示 t=0时 y的取值 即初值 %h表示步长 %f 表示函数值 %输出y表示方程的响应y y=y0+h*f2; end自定义函数f1.m function f= f1(x,y) f = -y+x+1;%微分方程 dy/dt=-y+t+1 初值y(0)=1 微分方程右边的剩余部分构成的函数 end运行结果  Euler法误差平方:
wucha_euler =
列 1 至 12
0 0.0110 0.0097 0.0086 0.0076 0.0067 0.0058 0.0051 0.0044 0.0039 0.0034 0.0029
列 13 至 21
0.0025 0.0022 0.0019 0.0016 0.0014 0.0012 0.0010 0.0009 0.0007改进程序：修改步长h，h分别取值0.1 0.05 0.01 0.001，取值t为0到1，对比结果。 主程序 clc; clear all; close all; h = [0.1 0.05 0.01 0.001];%步长 for j = 1:length(h) y0 = 1;%初值 t = 0:h(j):1;%x范围 y = exp(-t)+t;%真解 n = length(t); numy = zeros(1,n); f2 = -y0; numy(1) = y0; %欧拉法计算 for i=2:n numy(i) = euler1(y0,h(j),f2); y0 = numy(i); f2 = f1(t(i),y0); end %每次因为步长不一样 所以不能用矩阵存结果 s{j,:} = numy;%引入元胞类型 Cell 能包含任何类型的数据，比如数值、字符串、逻辑值甚至是Cell自身。 T{j,:}=t; end %绘图 figure; plot(t,y,'r-','linewidth',1); hold on; plot(T{1,:},s{1,:},'b-','linewidth',1); plot(T{2,:},s{2,:},'g-','linewidth',1); plot(T{3,:},s{3,:},'k-','linewidth',1); plot(T{4,:},s{4,:},'m-','linewidth',1); xlabel('t'); grid on; title('Euler法求系统的输出响应'); ylabel('输出响应y(t)'); legend('真解','Euler法 h=0.1','Euler法 h=0.05','Euler法 h=0.0.01','Euler法 h=0.001','location','northwest');自定义函数f1.m function f= f1(x,y) f = -y+x+1;%微分方程 dy/dt=-y+t+1 初值y(0)=1 微分方程右边的剩余部分构成的函数 end自定义函数euler1.m function y = euler1(y0,h,f2) %%输入参数 y0表示 t=0时 y的取值 即初值 %h表示步长 %f 表示函数值 %输出y表示方程的响应y y=y0+h*f2; end运行结果  作 者 | 郭志龙编 辑 | 郭志龙校 对 | 郭志龙 本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。 |
【本文地址】
今日新闻 |
推荐新闻 |