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

news2024/10/2 3:29:44

       目录

💥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.coloradmin.cn/o/338039.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

加盟管理系统挑选法则,看完不怕被坑!

经营服装连锁店铺究竟有多难?小编已经不止一次听到身边的老板,抱怨加盟连锁店铺难以管理了,但同时呢,也听到了很多作为加盟商的老板,抱怨总部给的支持和管理不到位。服装加盟店铺管理,到底有哪些难点呢&…

BFS广度优先遍历——Acwing 844. 走迷宫

1.BFS简介我们可以将bfs当做一个成熟稳重的人,一个眼观六路耳听八方的人,他每次搜索都是一层层的搜索,从第一层扩散到最后一层,BFS可以用来解决最短路问题。2.基本思想从初始状态S开始,利用规则,生成所有可…

window11 安装node及配置环境变量

一、安装环境 本教程演示的环境: 系统:win 11 64位 node.js下载地址: http://nodejs.cn/ node.js版本:长期支持版本(本教程基于16.15.0) 点击选中图标下载到电脑本地即可。 二、安装步骤 1、双击安装包,一…

华为10年经验测试工程师,整理出来的python自动化测试实战

前言 全书共分11章,第一章是基础,了selenium家谱,各种组件之间的关系以及一些必备知识。第二章告诉如何开始用python IDLE写程序以及自动化测试环境的搭建。第三章是webdriver API,我花了相当多时间对原先的文档,冗余…

HTML5之HTML基础学习笔记

列表标签 列表的应用场景 场景:在网页中按照行展示关联性的内容,如:新闻列表、排行榜、账单等特点:按照行的方式,整齐显示内容种类:无序列表、有序列表、自定义列表 这是老师PPT上的内容, 列表…

day10_面向对象基础

今日内容 零、 复习昨日 一、面向对象的概念 二、面向对象编程 三、内存图 零、 复习昨日 见晨考题 每日一数组题 写一个方法 用于合并两个int类型的数组 合并法则如下 {1,2,5,8,9}{1,3,0}---->{1,2,5,8,9,1,3,0} package com.qf.array;import java.util.Arrays;/*** --- 天…

基于Java+SpringBoot+Vue+uniapp前后端分离图书阅读系统设计与实现

博主介绍:✌全网粉丝3W,全栈开发工程师,从事多年软件开发,在大厂呆过。持有软件中级、六级等证书。可提供微服务项目搭建、毕业项目实战、项目定制✌ 博主作品:《微服务实战》专栏是本人的实战经验总结,《S…

MySQL在Linux上的四种安装方式

目录 前言 一、仓库安装 二、本地安装 三、容器安装 四、源码安装 前言 博主的配置信息: Windows版本:Win10 VMware虚拟机版本:Vmware Workstation Pro 17 Linux版本:Red Hat Enterprise Linux 9.1 MySQL版本:My…

一篇文章搞懂Cookie

目录 1 什么是Cookie 2 创建Cookie 3 浏览器查看Cookie 3.1 浏览器查看Cookie的第一种方式 3.2 浏览器查看Cookie的第二种方式 4 获取Cookie 5 修改Cookie 6 Cookie编码与解码 6.1 创建带中文Cookie 6.2 读取带中文Cookie 6.3 获取中文Cookie请求效果 6.4 解决创建和…

grafana9 使用消息模板配置发送企业微信(wecom)

一、grafana9告警设置: 1、进入告警消息模板介面 2、grafana 消息模板设置 template name : API_msg_tpl #名字随便 {{ define "myalert" }} **警报时间:** {{ .StartsAt.Format "2006-01-02 15:04:05 " }} {{ if gt (len .Labels) 0 }}**…

毕业5年,从月薪3000到年薪40w,我掌握了那些核心技能?(建议收藏)

大家好,我是静静~~是一枚一线大厂的测试开发工程师很多读者私信问我,自己时间不短了,随着工作年限的不断增长,感觉自己的技术水平与自己的工作年限严重不符。想跳槽出去换个新环境吧,又感觉自己的能力达不到心仪公司的…

Python_pytorch

python_pytorch 小土堆pytotch学习视频链接 from的是一个个的包(package) import 的是一个个的py文件(file.py) 所使用的一般是文件中的类(.class) 第一步实例化所使用的类,然后调用类中的方法(def) Dataset 数据集处理 import os from PIL impo…

本地(window)使用alist和RaiDav网盘挂载

一、背景 百度网盘的限速可能会让你转战阿里云盘,但是阿里云盘的缺点在于不能分享,网络上的资源都是通过各类网盘来分享的,这样就会让你可能同时拥有不同网盘的账号。 那么我们有没有一款工具,可以将这些网盘资源聚合一下&#xf…

RMQ--区间最值问题(在更)

RMQ(Range Minimum/Maximum Query)RMQ解决的问题ST算法 O(nlogn)线段树例题数列区间最大值最敏捷的机器人天才的记忆Frequent values总结(ST和线段树对比)RMQ解决的问题 RMQ是一个解决多个区间最值查询的算法,即区间最值查询&…

MySQL 创建数据表

在创建数据库之后,接下来就要在数据库中创建数据表。所谓创建数据表,指的是在已经创建的数据库中建立新表。 创建数据表的过程是规定数据列的属性的过程,同时也是实施数据完整性(包括实体完整性、引用完整性和域完整性&#xff09…

LwIP系列--线程通信消息结构

一、目的如果有小伙伴移植过LwIP,那么你肯定知道在LwIP源码中tcp/ip协议栈是作为一个单独的线程运行的,那么就有这样一个问题,我们从mac外设上收到的以太网数据包是如何交给tcp/ip线程进行处理的,用户发送的数据又是如何经过协议栈…

不学Python迟早会被淘汰?Python真有这么好的前景?

最近几年Python编程语言在国内引起不小的轰动,有超越Java之势,本来在美国这个编程语言就是最火的,应用的非常非常的广泛,而Python的整体语言难度来讲又比Java简单的很多。尤其是在运维的应用中非常的广泛,所以之前出了…

Ubuntu20.04无线网卡驱动安装

文章目录一.未安装无线网卡驱动的Ubuntu20.04联网方式二.Ubuntu20.04无线网卡驱动安装UbuntuU盘启动盘安装好Ubuntu 20.04之后,发现没有无线网络,不过有线可以用。一.未安装无线网卡驱动的Ubuntu20.04联网方式 比较简单的就是直接拉一条网线进行连接&am…

【C语言】宏定义 结构体 枚举变量的用法

目录 一、数据类型 二、C语言宏定义 三、C语言typedef重命名 四、 #define与typedef的区别 五、结构体 六、枚举变量 补充学习一点STM32的必备基础知识 一、数据类型 二、C语言宏定义 关键字:#define 用途:用一个字符串代替一个数字,…

214 情人节来袭,电视剧 《点燃我温暖你》李峋同款 Python爱心表白代码,赶紧拿去用吧

大家好,我是徐公,六年大厂程序员经验,今天为大家带来的是动态心形代码,电视剧 《点燃我温暖你》同款的,大家赶紧看看,拿去向你心仪的对象表白吧,下面说一下灵感来源。 灵感来源 今天&#xff…