拉普拉斯方程的解matlab,急求用matlab编写解拉普拉斯方程的程序
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
这么处理狄拉克边界条件,我的精确解和数值解的误差很大部只哪错了
附程序
% examp1a2.m
clear all
clc
N=5;M=5;n=N*M;
L=2;h=2;
tic
c=sqrt((L/(N-1))^2+(h/(M-1))^2);
dm=2.5*c;
mu=0.25;
E=1e5;
beta=1*10^5;
for j=1:M
for i=1:N
x(i+(j-1)*N)=L/(N-1)*(i-1);
y(i+(j-1)*N)=h/(M-1)*(j-(M+1)/2);
end
end
t=[-0.8611363;-0.339880;0.339880;0.8611363];
G=[0.3478548;0.6521452;0.6521452;0.3478548];
for j=1:(M-1)
for i=1:(N-1)
for ky=1:4
for kx=1:4
xg(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=(x(i+1+(j-1)*N)-x(i+(j-1)*N))/2*t(kx)+(x(i+1+(j-1)*N)+x(i+(j-1)*N))/2;
yg(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=(y(i+j*N)-y(i+(j-1)*N))/2*t(ky)+(y(i+j*N)+y(i+(j-1)*N))/2;
Gx(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=G(kx)*(x(i+1+(j-1)*N)-x(i+(j-1)*N))/2;
Gy(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=G(ky)*(y(i+j*N)-y(i+(j-1)*N))/2;
end
end
end
end
%inputing P
for i=1:n;
P(i,1)=1;
P(i,2)=x(i);
P(i,3)=y(i);
P(i,4)=(x(i))^2;
P(i,5)=x(i)*y(i);
P(i,6)=(y(i))^2;
end
%calculating K
K=0;
for i=1:(4*4*(N-1)*(M-1))
for j=1:n
d(j)=sqrt((xg(i)-x(j))^2+(yg(i)-y(j))^2);
if d(j)<=dm
W(j,j)=(exp(-(d(j)/c)^2)-exp(-6.25))/(1-exp(-6.25));
Wx(j,j)=-exp(-(d(j)/c)^2)*2*(xg(i)-x(j))/c/c/(1-exp(-6.25));
Wy(j,j)=-exp(-(d(j)/c)^2)*2*(yg(i)-y(j))/c/c/(1-exp(-6.25));
else
W(j,j)=0;
Wx(j,j)=0;
Wy(j,j)=0;
end
end
p=[1;xg(i);yg(i);(xg(i))^2;xg(i)*yg(i);(yg(i))^2];
px=[0,1,0,2*xg(i),yg(i),0];
py=[0,0,1,0,xg(i),2*yg(i)];
A=P'*W*P;
Ax=P'*Wx*P;
Ay=P'*Wy*P;
B=P'*W;
Bx=P'*Wx;
By=P'*Wy;
ph=p'*inv(A)*B;
phx=px*inv(A)*B-p'*inv(A)*Ax*inv(A)*B+p'*inv(A)*Bx;
phy=py*inv(A)*B-p'*inv(A)*Ay*inv(A)*B+p'*inv(A)*By;
for j=1:n
phi(1,2*j-1)=ph(j);
phi(1,2*j)=0;
phi(2,2*j-1)=0;
phi(2,2*j)=ph(j);
Ca(1,2*j-1)=phx(j);
Ca(1,2*j)=0;
Ca(2,2*j-1)=0;
Ca(2,2*j)=phx(j);
Cb(1,2*j-1)=phy(j);
Cb(1,2*j)=0;
Cb(2,2*j-1)=0;
Cb(2,2*j)=phy(j);
end
K=K+Gx(i)*Gy(i)*(Ca'*Ca+Cb'*Cb);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:(M-1)
for ik=1:4
xgb(4*(j-1)+ik)=L;
ygb(4*(j-1)+ik)=(y((j+1)*N)-y(j*N))/2*t(ik)+(y((j+1)*N)+y(j*N))/2;
Gyb(4*(j-1)+ik)=G(ik)*(y((j+1)*N)-y(j*N))/2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:(4*(M-1))
for j=1:n
db(j)=sqrt((xgb(i)-x(j))^2+(ygb(i)-y(j))^2);
if db(j)<=dm
Wb(j,j)=(exp(-(db(j)/c)^2)-exp(-6.25))/(1-exp(-6.25));
else
Wb(j,j)=0;
end
end
pb=[1;xgb(i);ygb(i);(xgb(i))^2;xgb(i)*ygb(i);(ygb(i))^2];
Ab=P'*Wb*P;
Bb=P'*Wb;
phb=pb'*inv(Ab)*Bb;
for j=1:n
phib(1,2*j-1)=phb(j);
phib(1,2*j)=0;
phib(2,2*j-1)=0;
phib(2,2*j)=phb(j);
end
K=K+beta*Gyb(i)*phib'*phib;
拉普拉斯方程的解matlab,急求用matlab编写解拉普拉斯方程的程序相关推荐
- 外推法程序matlab,急求用MATLAB用龙格库塔和外推法解一阶微分方程
共回答了21个问题采纳率:81% f=inline('-y+x+1','x','y'); %微分方程的右边项 dx=0.05; %x方向步长 xleft=0; %区域的左边界 xright=10; % ...
- matlab求曲线极值程序,matlab函数求极值matlab函数求极值.ppt
matlab函数求极值matlab函数求极值 * * 函数的极值 1.一元函数的极值 函数命令:fminbnd 调用格式:[x,feval,exitflag,output]=fminbnd(fun,x ...
- matlab ode45三体问题,小白急求~~关于ode45不能解的问题
ode45本身是变步长的呀,怎么设置最小步长?我改用其它ode函数计算都会出现这个警告Warning: Failure at t=5.373174e-001. Unable to meet inte ...
- matlab中求立方根,MATLAB基础入门
MATLAB有许多使用方法,但最基本,也是入门时首先要掌握的是MATLAB命令窗口(Command Window)的使用方法. MATLAB命令窗口是用于输入数据,运行MATLAB函数和脚本,并显示结 ...
- matlab trapz求二重积分,matlab求积分(超详细,含int integral integral2/3 quad trapz
matlab求积分 matlab求积分函数工具: int 用法1: 格式: int(fun,x,a,b) 功能: 计算定积分 用法2: 格式: int(f,x) 功能: 计算不定积分 注: 使用int ...
- matlab求解拉普拉斯方程,急求用matlab编写解拉普拉斯方程的程序
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 这么处理狄拉克边界条件,我的精确解和数值解的误差很大部只哪错了 附程序 % examp1a2.m clear all clc N=5;M=5;n=N*M; ...
- matlab求拉普拉斯方程,急求用matlab编写解拉普拉斯方程的程序
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 这么处理狄拉克边界条件,我的精确解和数值解的误差很大部只哪错了 附程序 % examp1a2.m clear all clc N=5;M=5;n=N*M; ...
- 【数学建模暑期培训】Matlab之求代数方程的符号解和数值解
文章目录 符号运算 1.1 符号对象的创建 1.2 代数方程的符号解 代数方程的数值解 线性方程组的数值解 非线性方程的数值解 符号运算 符号运算又称计算机代数,通俗地讲就是用计算机推导数学公式,如对 ...
- 矩阵特征值的用matlab,[急求]谁可以用matlab帮我运行求矩阵特征值的命令???...
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 v = Columns 1 through 6 0.0529 0.0228 - 0.0573i 0.0228 + 0.0573i 0.0222 + 0.0 ...
- matlab如何求极点,matlab求极点和零点
让每个人平等地提升自我一.实验目的 1.利用 MATLAB 的 fdatool 观察传递函数 H(z)的零极点分布.幅度响应.脉冲响 应及阶跃响应,并观察零极点分布与系统稳定性...... 假设为一个 ...
最新文章
- Owncloud-X安装配置
- 数学建模第五节2020.5.8-17补
- 深入理解Kafka(3)-Consumer
- 设置TextView文字
- GPT-GNN:图神经网络的生成式预训练 KDD 2020
- 实验4	C++程序的结构(4学时)
- Smarty的入门使用
- 在Redis集群技术上,你不可错过的四大集成者
- WebAssembly 介绍
- 语音识别系统原理介绍---从gmm-hmm到dnn-hmm
- android9.0 从driver到APP(2)--hardware
- Book-Manager 图书管理系统(基于SpringBoot、MyBatis)
- 智能中医诊疗系统php代码,智能新型中西医处方系统
- windows如何去除桌面图标箭头
- kubectl源码分析之rollout undo
- URLDownloadToFile缓存问题
- JavaSE基础知识(附上代码实现)1
- 有关于投资最优化的模型求解
- abaqus inp扫盲与提高 *MATRIX GENERATE,STIFFNESS的验证
- chrome恐龙游戏java_自动玩Chrome浏览器的小恐龙游戏