查看: 712|回复: 20|关注: 0

[已解决] 在使用ode45时,二阶微分方程无法直接手打输入,是通过for循环得到的式子,该怎么表达方程呢?

[复制链接]

新手

36 麦片

财富积分


050


3

主题

15

帖子

0

最佳答案
在求微分方程的时候,我的公式是通过一系列运算得到的,所以没办法直接输入到 ode45里面或者dx里面 ,该怎么解决




                               
登录/注册后可看大图





function dx=odeode(t,x)

syms t  u(t)  x1 x2

Af1=basicF3(x1,x2,u);       %Af1为状态函数1的基函数
n1=size(Af1,2);           %n1为基函数的维度:列数
A1=rand(n1,1);            %随机生成A1
A2=rand(n1,1);            %随机生成A2
A=[A1', A2']';            %A:状态方程的系数矩阵

Bf=basicF2(x1,x2);       %Bf为控制函数的基函数
n2=size(Bf,2);          %n2为基函数的维度:列数
B=rand(n2,1);           %B为控制函数u的系数矩阵
BF=Bf;                  %参考AF
U=sum(BF*B);              %控制策略方程

Afu1=basicF3(x1,x2,U) ;      %Af1为状态函数的基函数(代入u)
Afu2=Afu1;
AFU=[Afu1 zeros(1,size(Afu2,2));zeros(1,size(Afu1,2)) Afu2] ;
y=AFU*A
Y=[sum(y.*[1;0]);sum(y.*[0;1])]

%dx=@(t,x)[sum(y.*[1;0]);sum(y.*[0;1])]
dx=@(t,x)[subs(sum(y.*[1;0]),{x1,x2},{x(1),x(2)});subs(sum(y.*[0;1]),{x1,x2},{x(1),x(2)})]

end


                               
登录/注册后可看大图


function F=basicF(x1,x2,u)
%
%三元基函数
n1=1;
n2=1;
n3=1;
i=1;
for v3=0:n3
    for v2=0:n2
        for v1=0:n1
            a=nchoosek(n1,v1);
            b=nchoosek(n2,v2);
            c=nchoosek(n3,v3);
            p1=a*x1^v1*(1-x1)^(n1-v1);
            p2=b*x2^v2*(1-x2)^(n2-v2);
            p3=c*u^v3*(1-u)^(n3-v3);
            p=p1*p2*p3;
            F(i)=p;
            i=i+1;
            %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)*c*x3^v3*(1-x3)^(n3-v3)
        end
    end
end

                               
登录/注册后可看大图


function F=basicF2(x1,x2)
%
%二元基函数
n1=1;
n2=1;
i=1;
for v2=0:n2
    for v1=0:n1
        a=nchoosek(n1,v1);
        b=nchoosek(n2,v2);
        p1=a*x1^v1*(1-x1)^(n1-v1);
        p2=b*x2^v2*(1-x2)^(n2-v2);
        p=p1*p2;
        F(i)=p;
        i=i+1;
        %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)
    end
end

                               
登录/注册后可看大图



                               
登录/注册后可看大图


clear all
clc
tspan = [0,15];
x0 = [0;1];
[t, x] = ode45(@odeode, tspan, x0);

plot(t,x(:,1),'-',t,x(:,2),'-.')
legend('x1','x2')


                               
登录/注册后可看大图

出现的问题

                               
登录/注册后可看大图




回复主题 已获打赏: 0 积分

举报

论坛优秀回答者

10

主题

1698

帖子

370

最佳答案
  • 关注者: 84
发表于 2020-8-31 16:11:20 | 显示全部楼层
你好,这个问题我之前没有注意,只以为运算时间久一点,刚仔细看了一下,问题出在
  1. A1=rand(n1,1);            %随机生成A1
  2. A2=rand(n1,1);            %随机生成A2
复制代码

个人理解A1,A2随机生成后应在每次循环中保持不变,但在这里没有设置,所以每次会产生不同的迭代结果,导致无法收敛
代码简单修改如下
  1. clear all
  2. clc
  3. tspan = [0,15];
  4. x0 = [0;1];
  5. [t, x] = ode45(@odeode, tspan, x0);

  6. plot(t,x(:,1),'-',t,x(:,2),'-.')
  7. legend('x1','x2')

  8. function dx=odeode(t,x)

  9. syms u(t)  x1 x2

  10. Af1=basicF3(x1,x2,u);       %Af1为状态函数1的基函数
  11. n1=size(Af1,2);           %n1为基函数的维度:列数
  12. rng('default');
  13. A1=rand(n1,1);           %随机生成A1
  14. A2=rand(n1,1);            %随机生成A2
  15. A=[A1', A2']';            %A:状态方程的系数矩阵

  16. Bf=basicF2(x1,x2);       %Bf为控制函数的基函数
  17. n2=size(Bf,2);          %n2为基函数的维度:列数
  18. B=rand(n2,1);           %B为控制函数u的系数矩阵
  19. BF=Bf;                  %参考AF
  20. U=sum(BF*B);              %控制策略方程

  21. Afu1=basicF3(x1,x2,U) ;      %Af1为状态函数的基函数(代入u)
  22. Afu2=Afu1;
  23. AFU=[Afu1 zeros(1,size(Afu2,2));zeros(1,size(Afu1,2)) Afu2] ;
  24. y=AFU*A;
  25. Y=[sum(y.*[1;0]);sum(y.*[0;1])];

  26. %dx=@(t,x)[sum(y.*[1;0]);sum(y.*[0;1])]
  27. dx=[subs(sum(y.*[1;0]),{x1,x2},{x(1),x(2)});subs(sum(y.*[0;1]),{x1,x2},{x(1),x(2)})];
  28. dx=double(dx);
  29. end
  30. function F=basicF3(x1,x2,u)
  31. %
  32. %三元基函数
  33. n1=1;
  34. n2=1;
  35. n3=1;
  36. i=1;
  37. for v3=0:n3
  38.     for v2=0:n2
  39.         for v1=0:n1
  40.             a=nchoosek(n1,v1);
  41.             b=nchoosek(n2,v2);
  42.             c=nchoosek(n3,v3);
  43.             p1=a*x1^v1*(1-x1)^(n1-v1);
  44.             p2=b*x2^v2*(1-x2)^(n2-v2);
  45.             p3=c*u^v3*(1-u)^(n3-v3);
  46.             p=p1*p2*p3;
  47.             F(i)=p;
  48.             i=i+1;
  49.             %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)*c*x3^v3*(1-x3)^(n3-v3)
  50.         end
  51.     end
  52. end
  53. end
  54. function F=basicF2(x1,x2)
  55. %
  56. %二元基函数
  57. n1=1;
  58. n2=1;
  59. i=1;
  60. for v2=0:n2
  61.     for v1=0:n1
  62.         a=nchoosek(n1,v1);
  63.         b=nchoosek(n2,v2);
  64.         p1=a*x1^v1*(1-x1)^(n1-v1);
  65.         p2=b*x2^v2*(1-x2)^(n2-v2);
  66.         p=p1*p2;
  67.         F(i)=p;
  68.         i=i+1;
  69.         %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)
  70.     end
  71. end
  72. end
复制代码
1.PNG
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

10

主题

1698

帖子

370

最佳答案
  • 关注者: 84
发表于 2020-8-30 20:03:01 | 显示全部楼层
代码不全,比如basicF3函数未定义
回复此楼 已获打赏: 0 积分

举报

新手

36 麦片

财富积分


050


3

主题

15

帖子

0

最佳答案
 楼主| 发表于 2020-8-30 20:46:30 | 显示全部楼层
20141303 发表于 2020-8-30 20:03
代码不全,比如basicF3函数未定义

就是这一段function F=basicF(x1,x2,u)
少打了一个3,function F=basicF3(x1,x2,u), 用这个替换一下就是了,感谢
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

10

主题

1698

帖子

370

最佳答案
  • 关注者: 84
发表于 2020-8-30 22:04:02 | 显示全部楼层 |此回复为最佳答案
仅供参考,运行时间太长,我就不算了
  1. clear all
  2. clc
  3. tspan = [0,15];
  4. x0 = [0;1];
  5. [t, x] = ode45(@odeode, tspan, x0);

  6. plot(t,x(:,1),'-',t,x(:,2),'-.')
  7. legend('x1','x2')

  8. function dx=odeode(t,x)

  9. syms u(t)  x1 x2

  10. Af1=basicF3(x1,x2,u);       %Af1为状态函数1的基函数
  11. n1=size(Af1,2);           %n1为基函数的维度:列数
  12. A1=rand(n1,1);            %随机生成A1
  13. A2=rand(n1,1);            %随机生成A2
  14. A=[A1', A2']';            %A:状态方程的系数矩阵

  15. Bf=basicF2(x1,x2);       %Bf为控制函数的基函数
  16. n2=size(Bf,2);          %n2为基函数的维度:列数
  17. B=rand(n2,1);           %B为控制函数u的系数矩阵
  18. BF=Bf;                  %参考AF
  19. U=sum(BF*B);              %控制策略方程

  20. Afu1=basicF3(x1,x2,U) ;      %Af1为状态函数的基函数(代入u)
  21. Afu2=Afu1;
  22. AFU=[Afu1 zeros(1,size(Afu2,2));zeros(1,size(Afu1,2)) Afu2] ;
  23. y=AFU*A;
  24. Y=[sum(y.*[1;0]);sum(y.*[0;1])];

  25. %dx=@(t,x)[sum(y.*[1;0]);sum(y.*[0;1])]
  26. dx=[subs(sum(y.*[1;0]),{x1,x2},{x(1),x(2)});subs(sum(y.*[0;1]),{x1,x2},{x(1),x(2)})];
  27. dx=double(dx);
  28. end
  29. function F=basicF3(x1,x2,u)
  30. %
  31. %三元基函数
  32. n1=1;
  33. n2=1;
  34. n3=1;
  35. i=1;
  36. for v3=0:n3
  37.     for v2=0:n2
  38.         for v1=0:n1
  39.             a=nchoosek(n1,v1);
  40.             b=nchoosek(n2,v2);
  41.             c=nchoosek(n3,v3);
  42.             p1=a*x1^v1*(1-x1)^(n1-v1);
  43.             p2=b*x2^v2*(1-x2)^(n2-v2);
  44.             p3=c*u^v3*(1-u)^(n3-v3);
  45.             p=p1*p2*p3;
  46.             F(i)=p;
  47.             i=i+1;
  48.             %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)*c*x3^v3*(1-x3)^(n3-v3)
  49.         end
  50.     end
  51. end
  52. end
  53. function F=basicF2(x1,x2)
  54. %
  55. %二元基函数
  56. n1=1;
  57. n2=1;
  58. i=1;
  59. for v2=0:n2
  60.     for v1=0:n1
  61.         a=nchoosek(n1,v1);
  62.         b=nchoosek(n2,v2);
  63.         p1=a*x1^v1*(1-x1)^(n1-v1);
  64.         p2=b*x2^v2*(1-x2)^(n2-v2);
  65.         p=p1*p2;
  66.         F(i)=p;
  67.         i=i+1;
  68.         %p=a*x1^v1*(1-x1)^(n1-v1)*b*x2^v2*(1-x2)^(n2-v2)
  69.     end
  70. end
  71. end
复制代码
回复此楼 已获打赏: 2 积分

举报

新手

36 麦片

财富积分


050


3

主题

15

帖子

0

最佳答案
 楼主| 发表于 2020-8-30 22:38:43 | 显示全部楼层
20141303 发表于 2020-8-30 22:04
仅供参考,运行时间太长,我就不算了

感谢回答:handshake:loveliness:我运行了下
如果单独把调用的函数放在.m文件中,会直接报错
错误使用 double
无法从 function_handle 转换为 double。

不知道为什么放在一起运行就会  显示正忙
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

10

主题

1698

帖子

370

最佳答案
  • 关注者: 84
发表于 2020-8-31 08:50:07 | 显示全部楼层
将自定义函数另存.m文件,和放到脚本文件运行是一样的,都是会无报错运行,只是目前来看运行时间很长
回复此楼 已获打赏: 0 积分

举报

新手

36 麦片

财富积分


050


3

主题

15

帖子

0

最佳答案
 楼主| 发表于 2020-8-31 15:37:37 | 显示全部楼层
20141303 发表于 2020-8-31 08:50
将自定义函数另存.m文件,和放到脚本文件运行是一样的,都是会无报错运行,只是目前来看运行时间很长 ...

感谢,今早运行了好久一直没能出结果,在运行的时候我让y一直显示输出,发现y的方程系数是一直在变的,实际上用调用的basicF3和basicF2函数生成一个y之后,y应该是不变的,然后就y这个方程求微分。
然后我注释掉调用的函数,给y直接复制了一个之前运行时产生的一个方程,就出结果了,所以你的修改是没有报错的,感谢啊ღ( ´・ᴗ・` )比心

还有个问题想请教下你,为什么y的方程会一直变化呢?是不是就导致一直出不了结果?
回复此楼 已获打赏: 0 积分

举报

新手

36 麦片

财富积分


050


3

主题

15

帖子

0

最佳答案
 楼主| 发表于 2020-8-31 16:35:36 | 显示全部楼层
20141303 发表于 2020-8-31 16:11
你好,这个问题我之前没有注意,只以为运算时间久一点,刚仔细看了一下,问题出在

个人理解A1,A2随机生成 ...

ღ( ´・ᴗ・` )比心,千言万语汇成俩字,感谢!
还有个问题,得到的数值解x可以表达成方程吗,怎么能把结果代入到接下来的方程中?
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

10

主题

1698

帖子

370

最佳答案
  • 关注者: 84
发表于 2020-8-31 16:56:10 | 显示全部楼层
额,你是说微分方程的解析解?
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /5 下一条

快速回复 返回顶部 返回列表