展开全部
下面是我对这道题思路:
1:
对于第一题的微分方程组在matlab中输入[f,g]=dsolve('Df=exp(f*sin(t)+g)','Dg=exp(g*cos(t))+f','t')答案是 Explicit solution could not be found,即没有解析解,因此只能求其近似解,也就是数值解。
求解微分方程的数值解需要知道“函数的初值”“自变量的范围”,即f(0)=?,g(0)=?,变量t的取值范围是多少,然后迭代运算,得到在给定定义域内的近似值。
2:
差分法也就是我们知道的欧拉法(Euler)"思想是用先前的差商近似代替倒数",直白一些的编程说法即:f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长,若精度要求不高,一般可取0.01。在定义区间内迭代求解即可。
龙哥库塔法一般用于高精度的求解,即高阶精度的改进欧拉法。
常用的是四阶龙哥库塔:
编程语言如下:
y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;
k1=f(xi,yi)
k2=f(xi+h/2,yi+h*k1/2);
k3=f(xi+h/2,yi+h*k2/2);
k4=f(xi+h,yi+h*k3);
设置终止条件迭代求解。
思想就是这样,代码如下:
%% 龙哥库塔or欧拉法求解微分方程
t=0:0.01:3; %自变量范围
f = [];
g = [];
f(1) = 0.1; %f初值
g(1) = 1; %g初值
h=0.001;
for i=1:length(t)
%% 欧拉法
% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));
% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));
%% 龙哥库塔
kf1 = exp(f(i)*sin(t(i)))+g(i);
kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);
kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);
kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);
f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;
kg1 = exp(f(i)*cos(t(i)))+f(i);
kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);
kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);
kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);
g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;
end
figure(1)
plot(t,f(1:length(t)),t,g(1:length(t)));
legend('f函数','g函数')
1:
对于第一题的微分方程组在matlab中输入[f,g]=dsolve('Df=exp(f*sin(t)+g)','Dg=exp(g*cos(t))+f','t')答案是 Explicit solution could not be found,即没有解析解,因此只能求其近似解,也就是数值解。
求解微分方程的数值解需要知道“函数的初值”“自变量的范围”,即f(0)=?,g(0)=?,变量t的取值范围是多少,然后迭代运算,得到在给定定义域内的近似值。
2:
差分法也就是我们知道的欧拉法(Euler)"思想是用先前的差商近似代替倒数",直白一些的编程说法即:f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长,若精度要求不高,一般可取0.01。在定义区间内迭代求解即可。
龙哥库塔法一般用于高精度的求解,即高阶精度的改进欧拉法。
常用的是四阶龙哥库塔:
编程语言如下:
y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;
k1=f(xi,yi)
k2=f(xi+h/2,yi+h*k1/2);
k3=f(xi+h/2,yi+h*k2/2);
k4=f(xi+h,yi+h*k3);
设置终止条件迭代求解。
思想就是这样,代码如下:
%% 龙哥库塔or欧拉法求解微分方程
t=0:0.01:3; %自变量范围
f = [];
g = [];
f(1) = 0.1; %f初值
g(1) = 1; %g初值
h=0.001;
for i=1:length(t)
%% 欧拉法
% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));
% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));
%% 龙哥库塔
kf1 = exp(f(i)*sin(t(i)))+g(i);
kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);
kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);
kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);
f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;
kg1 = exp(f(i)*cos(t(i)))+f(i);
kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);
kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);
kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);
g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;
end
figure(1)
plot(t,f(1:length(t)),t,g(1:length(t)));
legend('f函数','g函数')
推荐律师服务:
若未解决您的问题,请您详细描述您的问题,通过百度律临进行免费专业咨询