当前位置： 首页 > news >正文

最小二乘支持向量机”在学习偏微分方程 （PDE） 解方面的应用（Matlab代码实现）

news 2025/6/26 10:33:47

       目录

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

💥1 概述

本代码说明了“最小二乘支持向量机”在学习偏微分方程 （PDE） 解方面的应用。提供了一个示例，并将获得的结果与精确的解决方案进行比较。

📚2 运行结果

主函数部分代码：

clc; clear all; close all

warning('off','all')

a0=0;

b0=1;

n=11;

h=(b0-a0)/n;

[X1,Y1]=meshgrid(a0+h:h:b0-h);

W=[];

for i=1:size(X1,2)

    Z=[X1(:,i),Y1(:,1)];

    W=[W ; Z];

end

subplot(2,3,1)

plot(W(:,1),W(:,2),'o')

hold on

[X,Y]=meshgrid(a0:h:b0);

W2=[];

for i=1:size(X,2)

    Z=[X(:,i),Y(:,1)];

    W2=[W2 ; Z];

end

L1=[];

for i=1:n+1

    L1=[L1 ; W2(i,:)];

end

L2=[];

for i=n*(n+1)+1:size(W2,1)

    L2=[L2 ; W2(i,:)];

end

L3=[L1(:,2) L1(:,1)];

L4=[L2(:,2) L2(:,1)];

plot(L1(:,1),L1(:,2),'s')

plot(L2(:,1),L2(:,2),'o')

plot(L3(:,1),L3(:,2),'p')

plot(L4(:,1),L4(:,2),'+')

title('Training points','Fontsize',14)

xlabel('x')

ylabel('y')

%% 

f=@(s,v) exp(-s).*(s-2+v.^3+6*v); % right hand side of the given PDE

gamma=10^14; % the regularization parameter

sig=0.95;  % kernel bandwidth

K=KernelMatrix(W,'RBF_kernel',sig);

x=W(:,1);

y=W(:,2);

xx1=x*ones(1,size(x,1));

xx2=x*ones(1,size(x,1));

cof1=2*(xx1-xx2')/(sig);

xx3=y*ones(1,size(y,1));

xx4=y*ones(1,size(y,1));

cof2=2*(xx3-xx4')/(sig);

Kxx=(-2/sig)*K + (cof1.^2) .* K;

Kyy=(-2/sig)*K + (cof2.^2) .* K;

Kx2x2=(   ( 12/(sig^2) - (12/sig)* (cof1.^2) +  (cof1.^4) ) .*K);

Ky2y2=(   ( 12/(sig^2) - (12/sig)* (cof2.^2) +  (cof2.^4) ) .*K);

Kx2y2=(   ( 4/(sig^2) - (2/sig)* (cof1.^2) - (2/sig)* (cof2.^2)  +  (cof1.^2).*(cof2.^2)  ) .*K);

Ky2x2=(   ( 4/(sig^2) - (2/sig)* (cof1.^2) - (2/sig)* (cof2.^2)  +  (cof1.^2).*(cof2.^2)  ) .*K);

K1T= Kx2x2+ Kx2y2 + Ky2x2+ Ky2y2;

m=size(K1T,1);

%*******************************************************************

KL1=KernelMatrix(W,'RBF_kernel',sig,L1);

L1b1x=L1(:,1)*ones(1,size(x,1));

L1b2x=x*ones(1,size(L1(:,1),1));

cofL1x=-2*(L1b1x'-L1b2x)/(sig);

L1b1y=L1(:,2)*ones(1,size(y,1));

L1b2y=y*ones(1,size(L1(:,2),1));

cofL1y=-2*(L1b1y'-L1b2y)/(sig);

KL1xx=(-2/sig)*KL1 + (cofL1x.^2) .* KL1;

KL1yy=(-2/sig)*KL1 + (cofL1y.^2) .* KL1;

KL1T= KL1xx+ KL1yy;

%*************************************************

KL2=KernelMatrix(W,'RBF_kernel',sig,L2);

L2b1x=L2(:,1)*ones(1,size(x,1));

L2b2x=x*ones(1,size(L2(:,1),1));

cofL2x=-2*(L2b1x'-L2b2x)/(sig);

L2b1y=L2(:,2)*ones(1,size(y,1));

L2b2y=y*ones(1,size(L2(:,2),1));

cofL2y=-2*(L2b1y'-L2b2y)/(sig);

KL2xx=(-2/sig)*KL2 + (cofL2x.^2) .* KL2;

KL2yy=(-2/sig)*KL2 + (cofL2y.^2) .* KL2;

KL2T= KL2xx+ KL2yy;

%*************************************************

KL3=KernelMatrix(W,'RBF_kernel',sig,L3);

L3b1x=L3(:,1)*ones(1,size(x,1));

L3b2x=x*ones(1,size(L3(:,1),1));

cofL3x=-2*(L3b1x'-L3b2x)/(sig);

L3b1y=L3(:,2)*ones(1,size(y,1));

L3b2y=y*ones(1,size(L3(:,2),1));

cofL3y=-2*(L3b1y'-L3b2y)/(sig);

KL3xx=(-2/sig)*KL3 + (cofL3x.^2) .* KL3;

KL3yy=(-2/sig)*KL3 + (cofL3y.^2) .* KL3;

KL3T= KL3xx+ KL3yy;

%*************************************************

KL4=KernelMatrix(W,'RBF_kernel',sig,L4);

L4b1x=L4(:,1)*ones(1,size(x,1));

L4b2x=x*ones(1,size(L4(:,1),1));

cofL4x=-2*(L4b1x'-L4b2x)/(sig);

L4b1y=L4(:,2)*ones(1,size(y,1));

L4b2y=y*ones(1,size(L4(:,2),1));

cofL4y=-2*(L4b1y'-L4b2y)/(sig);

KL4xx=(-2/sig)*KL4 + (cofL4x.^2) .* KL4;

KL4yy=(-2/sig)*KL4 + (cofL4y.^2) .* KL4;

KL4T= KL4xx+ KL4yy;

%*************************************************

KL1L1=KernelMatrix(L1,'RBF_kernel',sig,L1);

KL2L1=KernelMatrix(L2,'RBF_kernel',sig,L1);

KL3L1=KernelMatrix(L3,'RBF_kernel',sig,L1);

KL4L1=KernelMatrix(L4,'RBF_kernel',sig,L1);

%*************************************************

KL1L2=KernelMatrix(L1,'RBF_kernel',sig,L2);

KL2L2=KernelMatrix(L2,'RBF_kernel',sig,L2);

KL3L2=KernelMatrix(L3,'RBF_kernel',sig,L2);

KL4L2=KernelMatrix(L4,'RBF_kernel',sig,L2);

%************************************************

KL1L3=KernelMatrix(L1,'RBF_kernel',sig,L3);

KL2L3=KernelMatrix(L2,'RBF_kernel',sig,L3);

KL3L3=KernelMatrix(L3,'RBF_kernel',sig,L3);

KL4L3=KernelMatrix(L4,'RBF_kernel',sig,L3);

%************************************************

KL1L4=KernelMatrix(L1,'RBF_kernel',sig,L4);

KL2L4=KernelMatrix(L2,'RBF_kernel',sig,L4);

KL3L4=KernelMatrix(L3,'RBF_kernel',sig,L4);

KL4L4=KernelMatrix(L4,'RBF_kernel',sig,L4);

%************************************************

A= [K1T+1/gamma*eye(m) , KL1T , KL2T, KL3T , KL4T , zeros((n-1)^2,1) ;....

    KL1T' , KL1L1' , KL2L1' , KL3L1' , KL4L1' , ones(n+1,1) ;...

    KL2T' , KL1L2' , KL2L2' , KL3L2' , KL4L2' , ones(n+1,1) ;...

    KL3T' , KL1L3' , KL2L3' , KL3L3' , KL4L3' , ones(n+1,1) ;...

    KL4T' , KL1L4' , KL2L4' , KL3L4' , KL4L4' , ones(n+1,1) ;...

    zeros((n-1)^2,1)' , ones(n+1,1)' , ones(n+1,1)' , ones(n+1,1)' , ones(n+1,1)' , 0 ];

B=[f(W(:,1),W(:,2)); L1(:,2).^3 ; (1+L2(:,2).^3)*exp(-1)  ;  L3(:,1).*exp(-L3(:,1)) ; exp(-L4(:,1)).*(L4(:,1)+1) ; 0 ];

result=A\B;

alpha=result(1:m);

beta1=result(m+1:m+n+1);

beta2=result(m+n+2:m+2*n+2);

beta3=result(m+2*n+3:m+3*n+3);

beta4=result(m+3*n+4:m+4*n+4);

b=result(end);

%% Result for training points

yhat= (Kxx' + Kyy')* alpha + KL1 * beta1 + KL2* beta2 + KL3* beta3 + KL4* beta4 +b;

yexa=@(p,q) exp(-p).*(p+q.^3);

yexact=yexa(W(:,1),W(:,2));

Error1= yexact- yhat;

MAX_Absolute_error_training=max(abs(yhat-yexact));

RMSE_training=sqrt(mse(yhat-yexact));

fprintf('-------  training set ------------------\n\n')

fprintf('Max Abs Error on training set=%d\n',MAX_Absolute_error_training)

fprintf('RMSE on training set=%d\n\n',RMSE_training)

subplot(2,3,2)

plot3(W(:,1),W(:,2),yhat,'pr')

hold all

plot3(W(:,1),W(:,2),yexact,'sb')

title('Approximate and exact solution for training points','Fontsize',14)

xlabel('x')

ylabel('y')

zlabel('u')

NError=reshape(Error1,size(X1,1),size(Y1,1));

Xn=linspace(0,1,n-1);

Yn=linspace(0,1,n-1);

subplot(2,3,3)

surface(Xn,Yn,NError)

shading interp

xlabel('y','Fontsize',14)

ylabel('x','Fontsize',14)

set(gca,'Fontsize',20)

grid on

h=colorbar;

set(h,'fontsize',14);

title('Absolute errors for training set','Fontsize',14)

%% Result for test points

a0=0;

b0=1;

n=31;

h=(b0-a0)/n;

[X2,Y2]=meshgrid(a0+h:h:b0-h);

WT=[];

for i=1:size(X2,2)

    Z=[X2(:,i),Y2(:,1)];

    WT=[WT ; Z];

end

subplot(2,3,4)

plot(WT(:,1),WT(:,2),'o')

title('Test points','Fontsize',14)

xlabel('x')

ylabel('y')

Kt=KernelMatrix(W,'RBF_kernel',sig,WT);

xt=WT(:,1);

yt=WT(:,2);

xx1t=x*ones(1,size(xt,1));

xx2t=xt*ones(1,size(x,1));

cof1t=-2*(xx1t-xx2t')/(sig);

xx3t=y*ones(1,size(yt,1));

xx4t=yt*ones(1,size(y,1));

cof2t=-2*(xx3t-xx4t')/(sig);

Ktestxx=(-2/sig)*Kt + (cof1t.^2) .* Kt;

Ktestyy=(-2/sig)*Kt + (cof2t.^2) .* Kt;

KKlte1=KernelMatrix(WT,'RBF_kernel',sig,L1);

KKlte2=KernelMatrix(WT,'RBF_kernel',sig,L2);

KKlte3=KernelMatrix(WT,'RBF_kernel',sig,L3);

KKlte4=KernelMatrix(WT,'RBF_kernel',sig,L4);

Ytest= (Ktestxx' + Ktestyy')* alpha + KKlte1 * beta1 + KKlte2* beta2 + KKlte3* beta3 + KKlte4* beta4 + b;

yextest=yexa(WT(:,1),WT(:,2));

subplot(2,3,5)

plot3(WT(:,1),WT(:,2),Ytest,'pr')

hold on

plot3(WT(:,1),WT(:,2),yextest,'sb')

title('Approximate and exact solution for test points','Fontsize',14)

xlabel('x')

ylabel('y')

zlabel('u')

yextest=yexa(WT(:,1),WT(:,2));

MAX_Absolute_error_test=max(abs(Ytest-yextest));

RMSE_test=sqrt(mse(Ytest-yextest));

fprintf('-------  test set ------------------\n\n')

fprintf('Max Abs Error on test set=%d\n',MAX_Absolute_error_test)

fprintf('RMSE on test set=%d\n\n',RMSE_test)

fprintf('-------  Finished -----------------------\n\n')

Error= Ytest - yextest ;

Ytnew=reshape(Ytest,size(X2,1),size(Y2,1));

Ytexa=reshape(yextest,size(X2,1),size(Y2,1));

NError=reshape(Error,size(X2,1),size(Y2,1));

Xn=linspace(0,1,n-1);

Yn=linspace(0,1,n-1);

subplot(2,3,6)

surface(Xn,Yn,NError)

shading interp

xlabel('y','Fontsize',14)

ylabel('x','Fontsize',14)

set(gca,'Fontsize',20)

grid on

h=colorbar;

set(h,'fontsize',14);

title('Absolute errors for test set','Fontsize',14)

🎉3 参考文献

[1] Mehrkanoon S., Falck T., Suykens J.A.K., "Approximate Solutions to Ordinary Differential Equations Using Least Squares Support Vector Machines",IEEE Transactions on Neural Networks and Learning Systems, vol. 23, no. 9, Sep. 2012, pp. 1356-1367.

[2] Mehrkanoon S., Suykens J.A.K.,"LS-SVM approximate solution to linear time varying descriptor systems", Automatica, vol. 48, no. 10, Oct. 2012, pp. 2502-2511.

[3] Mehrkanoon S., Suykens J.A.K., "Learning Solutions to Partial Differential Equations using LS-SVM",Neurocomputing, vol. 159, Mar. 2015, pp. 105-116.

👨‍💻4 Matlab代码 

http://www.lryc.cn/news/2534.html

相关文章：

  • ISYSTEM调试实践8-winIDEA Analyzer功能1
  • 每日学术速递2.11
  • 宝塔搭建实战php开源likeadmin通用管理admin端vue3源码（二）
  • 网络基础-虚拟化工具-网桥
  • 剑指 Offer 14- II. 剪绳子 II
  • English Learning - Day55 作业打卡 2023.2.9 周四
  • pixhawk2.4.8-地面站配置-APM固件
  • golang 通道类型
  • 并发、并行、吞吐量、延迟、响应时间 含义理解
  • HTTP 和 HTTPS 的区别
  • 微搭低代码从入门到精通07-基础布局组件
  • Docker镜像的创建
  • 电子技术——MOS差分输入对
  • 树莓派 - 小记
  • 【论文解读|KDD2020】AKT. Context-Aware Attentive Knowledge Tracing
  • Geek Uninstaller：向流氓软件火力全开，超良心的软件彻底卸载工具
  • Java线程池
  • 2023-02-10 - 5 文本搜索
  • 华为OD机试 - 最近的医院（Python），简单直白
  • Leetcode.1223 掷骰子模拟
  • 数据分析到底该怎么学呢？讲真，真不难！
  • 活动星投票紫砂新青年制作一个投票活动
  • Git | 在IDEA中使用Git
  • ＜ Linux ＞：Linux 进程概念 (4)
  • 七、Java框架之MyBatisPlus
  • C语言柔性数组
  • 支付功能测试用例
  • 牛客网Python篇数据分析习题（一）
  • 【C语言】“指针类型”与“野指针”
  • Linux：软链接和硬链接的理解