matlab 怎样用龙格库塔法求二阶微分方程_百度知道
提问者采纳
初值要确定&假设a=b=c=d=1;&结果为：
非常感谢，但是a，b,c,d不是具体的值，怎么样求解析解？
二阶常系数其次线性微分方程y=dsolve('D2y+2/b*Dy*a^2*y','Dy(0)=c','y(0)=d') y = (b^(1/2)*tanh(a*b^(1/2)*((a^2*d^2 + b*c)/b)^(1/2)*(t/b + atanh((a*d)/(b^(1/2)*((a^2*d^2 + b*c)/b)^(1/2)))/(a*b^(1/2)*((a^2*d^2 + b*c)/b)^(1/2))))*((a^2*d^2 + b*c)/b)^(1/2))/a
怎样用四阶龙格库塔法求呢
四阶龙格库塔法是求数值解的
提问者评价
其他类似问题
龙格库塔法的相关知识
等待您来回答
下载知道APP
随时随地咨询
出门在外也不愁MATLABode45方法求解二阶微分方程 ，0.05(dx^2)/(d^2 t)+200x+0.8dx/dt-960+40288*P=0&br/&P的初始值为0.8。针对该微分方程采用数值解法，使用MATLAB的ode45（四阶五级的龙格库塔法）求解该常微分方程，并应用plot函数对解
MATLABode45方法求解二阶微分方程 ，0.05(dx^2)/(d^2 t)+200x+0.8dx/dt-960+40288*P=0P的初始值为0.8。针对该微分方程采用数值解法，使用MATLAB的ode45（四阶五级的龙格库塔法）求解该常微分方程，并应用plot函数对解 10
补充：并应用plot函数对解得的数值绘制其二维数据曲线图（P随t的变化）。求大神编程，谢谢。
不区分大小写匿名
相关知识等待您来回答
编程领域专家查看: 2442|回复: 5
龙格库塔法求解一阶微分方程组的程序
该用户从未签到
方程是4个一阶微分方程组.用龙格库塔法实现.最后画出Z1,Z2,Q1,Q2随时间演化的过程.
Z1'(t,Z1,Z2,Q1,Q2)=-4 / 3 * sqrt((1 + Z2 - Z1) * (1 + Z2 + 2 * Z1)) * sin(Q1) + 2 / 3 * sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) * sin(Q2)
Z2'(t,Z1,Z2,Q1,Q2)=2 / 3 * sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) * sin(Q1) - 4 / 3 * sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) * sin(Q2)
Q1'(t,Z1,Z2,Q1,Q2)=Z1 + (sqrt((1 + Z2 - Z1) * (1 + Z2 + 2 * Z1)) / (1 + Z2 - Z1) - sqrt((1 + Z2 - Z1) * (1 + Z2 + 2 * Z1)) / (1 + Z2 + 2 * Z1)) * cos(Q1) + sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) / (1 + Z2 - Z1) * cos(Q2)
Q2'(t,Z1,Z2,Q1,Q2)=Z2 + sqrt((1 + Z2 - Z1) * (1 + Z2 + 2 * Z1)) / (1 + Z2 - Z1) * cos(Q1) + (sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) / (1 - Z1 - 2 * Z2) - sqrt((1 + Z2 - Z1) * (1 - Z1 - 2 * Z2)) / (1 + Z2 - Z1)) *cos(Q2)
sqrt表示开方,
那位高手,帮小弟一把.
该用户从未签到
这个要自己慢慢研究了呵呵，楼主我也要做26个多的方程组成的方程组，也要用龙格库塔法，有机会可以多多交流下。要是你成功了，要拿出来分享下哈。
该用户从未签到
原帖由 chinamiracle 于
23:25 发表
这个要自己慢慢研究了呵呵，楼主我也要做26个多的方程组成的方程组，也要用龙格库塔法，有机会可以多多交流下。要是你成功了，要拿出来分享下哈。 那当然。
该用户从未签到
function main%主函数
tspan=[0:0.01:2];%求解时间为2秒也可以 tspan=[0 2]表示自动时间步长
x0=[0;0;0;0];%初值
[t,x]=ode45(@ff,tspan,x0);%用龙格库塔法求解
plot(t,x(:,1))
plot(t,x(:,2))
plot(t,x(:,3))
plot(t,x(:,4))
function dx=ff(t,x)
dx=zeros(4,1);
dx(1)=-4 / 3 * sqrt((1 + x(2) - x(1)) * (1 + x(2) + 2 * x(1))) * sin(x(3)) + 2 / 3 * sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) * sin(x(4));
dx(2)=2 / 3 * sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) * sin(x(3)) - 4 / 3 * sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) * sin(x(4));
dx(3)=x(1) + (sqrt((1 + x(2) - x(1)) * (1 + x(2) + 2 * x(1))) / (1 + x(2) - x(1)) - sqrt((1 + x(2) - x(1)) * (1 + x(2) + 2 * x(1))) / (1 + x(2) + 2 * x(1))) * cos(x(3)) + sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) / (1 + x(2) - x(1)) * cos(x(4));
dx(4)=x(2) + sqrt((1 + x(2) - x(1)) * (1 + x(2) + 2 * x(1))) / (1 + x(2) - x(1)) * cos(x(3)) + (sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) / (1 - x(1) - 2 * x(2)) - sqrt((1 + x(2) - x(1)) * (1 - x(1) - 2 * x(2))) / (1 + x(2) - x(1))) *cos(x(4));
随便帮楼主写了一个，还没装matlab没法试，应该不会有什么大问题
该用户从未签到
原帖由 zhuchangliang 于
10:05 发表
那当然。 这个问题其实我跟你发帖那天就解决了哈，只是那天没付好值结果老是出不来。
该用户从未签到
我大致成功了
f1in=[0 0.2]
f2in=[0 0.2]
f3in=[0 0]
f4in=[0 0]
f1=inline('-4/3*sqrt((1+Z2-Z1)*(1+Z2+2*Z1))*sin(Q1)+2/3*sqrt((1+Z2-Z1)*(1-Z1-2*Z2))*sin(Q2)','Z1','Z2', 'Q1','Q2')
f2=inline('2/3*sqrt((1+Z2-Z1)*(1+Z2+2*Z1))*sin(Q1)-4/3*sqrt((1+Z2-Z1)*(1-Z1-2*Z2))*sin(Q2)','Z1','Z2', 'Q1','Q2')
f3=inline('Z1+(sqrt((1+Z2-Z1)*(1+Z2+2*Z1))/(1+Z2-Z1)-sqrt((1+Z2-Z1)*(1+Z2+2*Z1))/(1+Z2+2*Z1))*cos(Q1)+sqrt((1+Z2-Z1)*(1-Z1-2*Z2))/(1+Z2-Z1)*cos(Q2)','Z1','Z2', 'Q1','Q2')
f4=inline('Z2-sqrt((1+Z2-Z1)*(1+Z2+2*Z1))/(1+Z2-Z1)*cos(Q1)+(sqrt((1+Z2-Z1)*(1-Z1-2*Z2))/(1-Z1-2*Z2)-sqrt((1+Z2-Z1)*(1-Z1-2*Z2))/(1+Z2-Z1))*cos(Q2)','Z1','Z2', 'Q1','Q2')
resf1(1,[1 2])=f1in
resf2(1,[1 2])=f2in
resf3(1,[1 2])=f3in
resf4(1,[1 2])=f4in
for i=1:iter
kf11=feval(f1, f1in(2),f2in(2),f3in(2),f4in(2))
kf21=feval(f2, f1in(2),f2in(2),f3in(2),f4in(2))
kf31=feval(f3, f1in(2),f2in(2),f3in(2),f4in(2))
kf41=feval(f4, f1in(2),f2in(2),f3in(2),f4in(2))
kf12=feval(f1,f1in(2)+0.5*h*kf11,f2in(2)+0.5*h*kf21, f3in(2)+0.5*h*kf31,f4in(2)+0.5*h*kf41)
kf22=feval(f2,f1in(2)+0.5*h*kf11,f2in(2)+0.5*h*kf21, f3in(2)+0.5*h*kf31,f4in(2)+0.5*h*kf41)
kf32=feval(f3,f1in(2)+0.5*h*kf11,f2in(2)+0.5*h*kf21, f3in(2)+0.5*h*kf31,f4in(2)+0.5*h*kf41)
kf42=feval(f4,f1in(2)+0.5*h*kf11,f2in(2)+0.5*h*kf21, f3in(2)+0.5*h*kf31,f4in(2)+0.5*h*kf41)
kf13=feval(f1,f1in(2)+0.5*h*kf12,f2in(2)+0.5*h*kf22, f3in(2)+0.5*h*kf32,f4in(2)+0.5*h*kf42)
kf23=feval(f2,f1in(2)+0.5*h*kf12,f2in(2)+0.5*h*kf22, f3in(2)+0.5*h*kf32,f4in(2)+0.5*h*kf42)
kf33=feval(f3,f1in(2)+0.5*h*kf12,f2in(2)+0.5*h*kf22, f3in(2)+0.5*h*kf32,f4in(2)+0.5*h*kf42)
kf43=feval(f4,f1in(2)+0.5*h*kf12,f2in(2)+0.5*h*kf22, f3in(2)+0.5*h*kf32,f4in(2)+0.5*h*kf42)
kf14=feval(f1,f1in(2)+h*kf13,f2in(2)+h*kf23, f3in(2)+h*kf33,f4in(2)+h*kf43)
kf24=feval(f2,f1in(2)+h*kf13,f2in(2)+h*kf23, f3in(2)+h*kf33,f4in(2)+h*kf43)
kf34=feval(f3,f1in(2)+h*kf13,f2in(2)+h*kf23, f3in(2)+h*kf33,f4in(2)+h*kf43)
kf44=feval(f4,f1in(2)+h*kf13,f2in(2)+h*kf23, f3in(2)+h*kf33,f4in(2)+h*kf43)
kkf1=f1in(2)+h*(kf11+2*kf12+2*kf13+kf14)/6
kkf2=f2in(2)+h*(kf21+2*kf22+2*kf23+kf24)/6
kkf3=f3in(2)+h*(kf31+2*kf32+2*kf33+kf34)/6
kkf4=f4in(2)+h*(kf41+2*kf42+2*kf43+kf44)/6
f1in(2)=kkf1
f2in(2)=kkf2
f3in(2)=kkf3
f4in(2)=kkf4
xx=f1in(1)+i*h
resf1(i+1,[1 2])=[xx kkf1]
resf2(i+1,[1 2])=[xx kkf2]
resf3(i+1,[1 2])=[xx kkf3]
resf4(i+1,[1 2])=[xx kkf4]
plot(resf1(:,1),resf1(:,2),'b',resf2(:,1),resf2(:,2),'R')
plot(resf3(:,1),resf3(:,2),'b',resf2(:,1),resf2(:,2),'R')
工作时间：8:00-24:00
百思(Baisi.net)急！求大神用Matlab四阶龙格库塔解个方程！ - 数学 - 小木虫 - 学术 科研 第一站
&& 查看话题
急！求大神用Matlab四阶龙格库塔解个方程！
万分感谢，求哪位大神用matlab四阶龙格库塔帮我解下面这个方程
dy/dx=(-1.588)*10^(-7)*y^(2/7)*(9.617y^(2/7)-5)^(1/2), y(0)=0.11
求看见的大神帮我解一下，拜托了，万分感谢！！！！！！
function SolveRK4
y0 = 0.11;
=ode45(@myfun,,y0);
plot(t,y,'linewidth',3);
xlabel('x');
ylabel('y(x)')
function dt = myfun(t,y)
untitled.jpg : Originally posted by laosam280 at
function SolveRK4
y0 = 0.11;
=ode45(@myfun,,y0);
plot(t,y,'linewidth',3);
xlabel('x');
ylabel('y(x)')
function dt = myfun(t,y)
untitled.jpg
... 感谢大神帮助，万分感谢
但是我的初始值搞错了，还有我想知道几个特定y值下对应的
麻烦大神给我重新算一下，拜托了
dy/dx=(-1.588)*10^(-7)*y^(2/7)*(9.617y^(2/7)-5)^(1/2), y(0)=0.1918
求y=0.096,0.093,0.09,0.087,0.079,0.072,0.065,0.056,0.052
拜托了大神！！！ : Originally posted by laosam280 at
function SolveRK4
y0 = 0.11;
=ode45(@myfun,,y0);
plot(t,y,'linewidth',3);
xlabel('x');
ylabel('y(x)')
function dt = myfun(t,y)
untitled.jpg
... 感谢大神，万分感谢
但是我初始值错了，还有就是，我想知道几个特定y值下对应的X值
dy/dx=(-1.588)*10^(-7)*y^(2/7)*(9.617y^(2/7)-5)^(1/2), y(0)=0.1918
求y=0.1,0.096,0.093,0.09,0.087，0.079,0.072,0.065,0.056,0.05,0.04,0.035时对应的X值，麻烦大神重新帮我算一下，拜托了，我不懂啊 你这个方程右端因为有一个很小的系数，使得方程的解在很长的时间内都稳定在初值附近。你的初值给的是0.1918，要想解减少到哪怕0.09，x会上好几十万这么大。我刚刚试了一下。
代码都给你了，你自己去matlab上运行就是了。检查一下你的方程，看看有没有什么错误。 : Originally posted by laosam280 at
你这个方程右端因为有一个很小的系数，使得方程的解在很长的时间内都稳定在初值附近。你的初值给的是0.1918，要想解减少到哪怕0.09，x会上好几十万这么大。我刚刚试了一下。
代码都给你了，你自己去matlab上运行就 ... haode ,xiexie
var cpro_id = 'u1216994';
欢迎监督和反馈：本帖内容由
提供，小木虫为个人免费站点，仅提供交流平台，不对该内容负责。欢迎协助我们监督管理，共同维护互联网健康，如果您对该内容有异议，请立即发邮件到
联系通知管理员，也可以通过QQ周知，我们的QQ号为：8835100
我们保证在1个工作日内给予处理和答复，谢谢您的监督。
小木虫，学术科研第一站，为中国学术科研研究提供免费动力
欢迎监督，发现不妥请立即
E-mail： & QQ：8835100下载此文档需消耗20财富值
您当前的财富值为2分，不足以支付下载此文档哦！
您可以通过以下途径获取财富值：
上传文档(更多财富值详情查看)
文件大小：KB所需财富值：1分
您当前拥有财富值20分
(在一个月内重复下载，不扣财富值)
您的收藏添加到了个人中心的
可能由于重复收藏或网络原因造成，请稍候重试
取消收藏成功
取消收藏失败
可能由于网络原因造成，请稍候重试
阅读排行榜
贡献者：606689abcdef
分类：理学
阅读次数：19334次
下载次数：1034次
下载所需财富值：2币
内容简介：
matlab编的4阶龙格库塔法解微分方程的程序
相关精品书推荐}